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Dear  Tom: 

Attached  to  this  letter  is  the  report  of  our  first  year's  efforts  under  Grant  No:  N00014-93-1- 
0763,  fflQtitled:  "Stochastic  Nonlinear  Dynamics  of  Floating  Structures".  This  work  was  carried  out  by 
me  and  ONR  Grant-supported  Graduate  Assistant,  Patrick  Bar-Avi.  In  addition,  related  work  was 
carried  out  by  Ed  Weinstein  and  me  on  the  issue  of  solving  for  the  moments  of  a  system  governed  by 
the  Fokker-Planck  equation. 

Attached  are  the  following  documents: 

1 .  Nonlinear  Dynamics  of  an  Articulated  Tower  in  the  Ocean,  submitted  for  review  and  publication 
to  the  Journal  of  Sound  and  Vibration 

2.  Nonlinear  Dynamics  of  an  Articulated  Tower  in  the  Ocean,  submitted  for  review,  presentation 
and  publication  at  the  ASME  Winter  Annual  Meeting 

3.  The  van  Kampen  Expansion  for  the  Fokker-Planck  Equation  of  a  Duffing  Oscillator,  accepted  for 
publication  in  Jmimal  of  Statistical  Physics 

4.  The  van  Kampen  Expansion  for  the  Fokker-Planck  Equation  of  a  Duffing  Oscillator  Excited  by 
Colored  Noise,  accepted  for  publication  in  Journal  of  Statistical  liiysics 

5.  The  van  Kampen  Expansion  for  a  Linked  Fokker-Planck  Equation  of  a  Diffing  Oscillator  Excited 
by  Colored  Noise,  submitted  for  review  and  publication  in  Journal  of  Statistical  niysics. 
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.  .*e  c‘'.sence  of  our  work  on  the  offshore  problem  is  embodied  in  document  1  above.  Our  approach  to 
the  problem  has  been,  and  continues  to  be,  based  on  an  examination  of  the  physics  of  the  environment 
tjid  its  interaction  with  the  structure.  This  allows  us  to  proceed  along  a  research  path  that  best 
examines  the  importance  of  each  force  component,  each  nonlinearity,  and  helps  us  decide  which 
terms  in  our  analytical  model  can  be  ignored  or  must  be  retained.  In  addition,  by  going  back  and  forth 
between  analytical  model  and  simulation,  we  are  able  to  estimate  the  loss  in  accuracy  resulting  from 
any  particular  approximation. 

It  is  in  this  manner  that  we  have  proceeded  with  the  work  you  see  before  you  on  the 
articulated  tower.  We  have  explored  various  behavior  regimes,  including  a  chaotic  one  and  a  friction 
damping  effect.  Our  main  conclusion  during  this  phase  is  that  a  simplified,  classical  nonlinear 
oscillator  model,  such  as  the  Duffing  or  the  van  dcr  Pol,  will  not  be  representative  of  the  highly 
complex  and  nonlinear  offshore  fluid-structure  interaction  problem.  Thus,  the  need  to  make  a 
connection  between  the  physical  world  and  the  mathematics  we  choose  to  model  it.  We  prefer  a  top- 
down  approach.  It  is  always  easier  to  a  posteriori  drop  terms  shown  to  be  negligible  than  to  add  terms 
required  to  better  fit  the  behavior. 

On  the  articulated  tower,  in  the  current,  second,  year  of  the  grant,  we  are  generalizing  the 
model  to  include  effects  such  as  current,  and  allowing  the  structure  to  move  freely  about  its  base 
support  hinge.  Thus,  transverse  forces  due  to  vortex  shedding  can  be  incorporated. 

Our  work  on  the  van  Kampen  expansion  of  the  Fokker-Planck  also  continues.  The  above  three 
manuscripts  (3,4,5)  focused  on  unforced  motion  about  an  equilibrium  due  to  noise.  This  permitted  us 
to  examine  the  method,  and  extend  it  to  include  colored  noise  models,  and  also  to  consider  coupled 
oscillator  models,  which  have  broad  applicability  in  the  physical  sciences  and  engineering.  Our 
coupled  model  was  of  a  Duffing  coupled  to  a  linear  harmonic  oscillator,  with  colored  noise  forcing. 
These  are  being  published  in  J.  Stat.  Phys.,  where  the  original  papers  by  van  Kampen  were  published. 
Our  next  step  is  to  consider  particular  forced  cases  for  the  above.  As  for  all  nonlinear  oscillations, 
each  loading  case  requires  specific  considerations. 

As  usual,  I  would  be  pleased  to  visit  you  to  provide  you  with  more  details,  and  you  are 
welcome  here  for  a  working  visit. 

Thank  you  for  your  interest  in  our  work.  Your  support  is  appreciated.  Best  regards. 


Sincerely  yours. 


Haym  Benaroya 

cc:  Administrative  Grants  Officer  QUALITY  INSPECTED  3 
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Abstract 


This  paper  presents  studies  on  the  response  of  an  articulated  tower  in  the  ocean  subjected  to  deterministic 
and  random  wave  loading.  The  tower  is  modeled  as  an  upright  rigid  pendulum  with  a  concentrated  mass  at 
the  top  and  having  one  angular  degree  of  freedom  about  a  hinge  with  coulomb  damping.  In  the  derivation  of 
the  differential  equation  of  motion,  nonlinear  terms  due  to  geometric  (large  angle)  and  fluid  forces  (drag  and 
inertia)  are  included.  The  wave  loading  is  derived  using  Morison’s  equation  in  which  the  velocity  and  acceleration 
of  the  fluid  are  determined  along  the  instantaneous  position  of  the  tower,  causing  the  equation  of  motion  to 
be  highly  nonlinear.  Furthermore,  since  the  differential  equation’s  coefficients  are  time-dependent  (periodic), 
parametric  instability  can  occur  dependiug  on  system  parameters  such  as  wave  height  and  frequency,  buoyancy, 
and  drag  coefficient.  The  nonlinear  differential  equation  is  then  solved  numericaUy  using  ‘ACSL’  software.  The 
response  of  the  tower  to  deterministic  wave  loading  is  investigated  and  a  stability  analysis  is  performed  (resonance 
and  parametric  instability).  To  solve  the  equation  for  random  loading,  the  Pierson-Moskowitz  power  spectrum, 
describing  the  wave  height,  is  first  transformed  into  an  approximate  time  history  using  Borgman’s  method  with 
slight  modification.  The  equation  of  motion  is  then  solved,  and  the  influence  on  the  tower  response  of  different 
parameter  values  such  as  buoyancy,  initial  conditions  and  wave  height  and  frequency,  is  investigated. 
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1  Review  and  Proble  )efinition 


Compliant  platforms  such  as  articulated  towers  are  economically  attractive  for  deep  water  conditions 
because  of  their  reduced  structural  weight  compared  to  conventional  platforms.  The  foundation  of 
the  tower  does  not  resist  lateral  forces  due  to  wind,  waves  and  currents;  instead,  restoring  moments 
are  provided  by  a  large  buoyancy  force,  a  set  of  guylines  or  a  combination  of  both.  These  structures 
have  a  fundamental  frequency  well  below  the  w®  ■  '■■‘r  frequency.  As  a  result  of  the  relatively 
large  displacements,  geometric  nonlinearity  is  an  impc  uit  consideration  in  the  analysis  of  such  a 
structure.  The  analysis  and  investigation  of  these  kind  of  problems  can  be  divided  into  two  major 
groups;  deterministic  and  random  wave  and/or  current  loading.  We  briefly  review  work  in  this  area 
in  the  next  two  subsections. 

1.1  Deterministic  loading 

Chakrabarti  and  Cotter  (1979)  [3]  analyzed  the  motion  of  articulated  tower.  The  tower  is  artic¬ 
ulated  by  a  imiversal  joint  having  single  degree  of  freedom.  They  assumed,  linear  waves,  small 
perturbations  about  an  equilibritun  position,  linear  drag  force  and  that  the  wind  current  and  wave 
are  colinear.  Their  resulting  equation  of  motion  is, 

-H  -I- CV’ =  (1) 

where  I  is  the  total  moment  of  inertia  including  added  mass,  is  the  nonlinear  drag  term.  Dip  is 
the  structural  damping,  Cipis  the  restoring  moment  due  to  buoyancy  and  Mo  is  the  wave  moment. 
An  analytical  solution  is  then  compared  to  experimental  results  to  show  good  agreement  as  long  as 
the  system  is  inertia  predominant,  and  not  drag  predominant. 

In  a  later  paper,  Chakrabarti  and  Cotter  (1980)  [4]  investigated  transverse  motion,  the  motion 


perpendicular  to  the  horizontal  velocity.  The  tower  pivot  is  assumed  to  have  two  angular  degrees 
of  freedom  and  is  taken  to  be  frictionless.  It  is  also  assumed  that  the  motion  is  not  coupled,  so 
the  inline  solution  is  obtained  (the  same  as  in  the  previous  paper)  from  which  the  relative  velocity 
between  the  tower  and  the  wave  is  obtained.  The  lift  force  (in  the  transverse  direction)  can  then 
be  obtained  and  the  linear  equation  of  motion  is  solved  analytically  and  compared  to  experimental 
results.  The  comparison  shows  good  agreement,  especially  when  the  drag  terms  are  small. 

Jain  and  Kirk  (1981)  [9]  investigated  the  dynamic  response  of  a  double  articulated  offshore 
structure  to  waves  and  current  loading.  They  assumed  four  degrees  of  freedom,  two  angular  for 
each  link.  The  equations  of  motion  were  derived  using  Lagrange  equations.  In  deriving  the  equations 
of  motion  the  following  assumptions  were  made:  drag  and  inertia  forces  tangent  to  the  tower  are 
negligible,  and  the  wave  and  current  velocities  are  evaluated  at  the  upright  position  (small  angles 
assumption).  The  linearized  equations  were  solved  to  find  the  natural  frequencies  of  the  system 
and  then  numerically  solved  to  find  the  response  due  to  colinear  and  non-colinear  current  and  wave 
velocities.  They  found  that  when  the  wave  and  the  current  velocities  are  colinear,  the  response  of 
the  top  is  sinusoidal,  while  for  noncolinear  velocities  the  response  is  a  complex  three  dimensional 
whirling  oscillation. 

Thompson  et  al.  (1984)  [16]  investigated  the  motions  of  an  articulated  mooring  tower.  They 
modeled  the  structure  2is  a  bilinear  oscillator  which  consists  of  two  linear  oscillators  having  different 
stiffnesses  for  each  half  cycle, 

mx  +  cx  +  (fci,  fc2)x  =  Fosinwf,  (2) 

where  ki,  k2  are  the  stiffnesses  for  x  >  0  and  x  <  0  respectively.  The  equation  is  solved  numerically 
for  different  spring  ratios  and,  as  expected,  harmonic  and  subharmonic  resonances  appeared  in  the 


3 


response.  A  comparison  between  the  response  and  experimental  results  of  a  reduced-scale  model 
showed  good  agreement  in  the  m2un  phenomenon. 

Choi  and  Lou  (1991)  [5]  have  investigated  the  behaviour  of  an  articulated  offshore  platform.  They 
modeled  it  as  an  upright  pendulum  having  one  degree  of  freedom,  with  linear  springs  at  the  top 
having  different  stiffnesses  for  positive  and  negative  displacements  (bilinear  oscillator).  The  equation 
of  motion  is  simplified  by  expanding  nonlinear  terms  into  a  power  series  and  retaining  only  the  first 
two  terms.  They  assumed  that  the  combined  drag  and  inertia  moment  is  a  harmonic  function.  The 
discontinuity  in  the  stiffness  is  assumed  to  be  small,  and  thus  replaced  by  an  equivalent  continuous 
function  using  a  least-squares  method  to  get  the  following  Duffing  equation 

16  +  c0  +  ki6  +  +  ks^  =  Mq  cos  uit,  (3) 

where  ibi,  Jfcj,  k^  are  spring  constants  depending  on  buoyancy,  gravity  and  the  mooring  lines.  The 
equation  of  motion  is  solved  analytically  and  numerically,  and  stability  analysis  is  performed.  The 
analytical  solution  agrees  very  well  with  the  numerical  solution.  The  main  results  of  their  analysis 
are  that  as  damping  decreases,  jump  phenomena  2uid  higher  subharmonics  occur,  and  chaotic  motion 
occurs  only  for  large  waves  and  near  the  first  subharmonic  (excitation  frequency  equals  twice  the 
fundamental  frequency);  the  system  is  very  sensitive  to  initial  conditions. 

Seller  and  Niedzwecki  (1992)  [14]  investigated  the  response  of  a  multi-articulated  tower  in  pl£in- 
ner  motion  (upright  multi-pendulum)  to  account  for  the  tower  flexibility.  The  restoring  moments 
(buoyancy  and  gravity)  were  taken  as  linear  rotational  springs  between  each  link,  although  the 
authors  say  that  nonlinear  springs  are  more  adequate  for  this  model.  Each  link  is  assumed  to  have 
a  different  cross  section  and  density.  The  equations  of  motion  are  derived  using  Lagrange’s  equa¬ 
tions,  in  which  the  generalized  coordinates  are  the  angular  deflections  of  each  link.  The  equations 
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I  in  matrix  form  are 

(M){«)  +  (Kl{9}  =  M,  (4) 

where  [M]  is  a  mass  matrix  that  includes  the  actual  mass  of  the  link  and  added  mass  terms,  while 

» 

the  stiffness  matrix  [/f]  includes  buoyancy  and  gravity  effects.  Damping  amd  drag  forces  are  not 
included  in  the  model.  The  homogeneous  equations  for  a  tri-articulated  tower  are  numerically 
•  solved  to  study  the  effects  of  different  parameters,  such  as  link  length,  material  density  and  spring 
stiffness,  on  the  natural  frequency  of  the  system. 

Gotlib  et  al.  (1992)  [7]  analyzed  the  nonlinear  response  of  a  single  degree  of  freedom  articulated 

I 

tower.  In  the  derivation  of  the  equation,  the  expressions  for  the  buoyancy  moment  arm,  added  mass 
term,  and  drag  and  inertia  moments  are  evaluated  along  the  stationary  upright  tower  position  and 
,  not  at  the  instantaneous  position  of  the  tower.  The  governing  equation  is  of  the  form 

0  +  7^  +  i2(^)  =  M(0,O,  (5) 

'  where  R{B)  =  asind  and  a  is  linear  function  of  buoyancy  and  gravity,  M{0,t)  is  the  drag  moment. 
Approximated  harmonic  and  subharmonic  solutions  are  derived  using  a  finite  Fourier  series  expan¬ 
sion,  and  stability  analysis  is  performed  by  a  Lyapunov  function  approach.  The  solution  shows  a 
jump  phenomenon  in  primary  2Uid  the  secondary  resonances. 

1.2  Random  loading 

Muhuri  and  Gupta  (1983)  [12]  investigated  the  stochastic  stability  of  a  buoyant  platform.  They 
used  a  linear  single  degree  of  freedom  model  as  follows 

X  -I-  2cx  -I-  (H-  G(f))x  =  0,  (6) 
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►  where  x  is  the  displacement,  c  is  the  damping  coefficient  and  G(t)  is  a  stochastic  time-dependent 
function  due  to  buoyancy.  It  is  assumed  that  6^(t)  is  a  narrow-band  random  process  with  zero-mean. 
A  criterion  for  the  mean  square  stability  is  obtained  from  which  the  following  results  are  found:  for 
c  >  1  the  system  is  always  stable,  and  for  0  <  c  <  1  there  are  regions  of  stability  and  instability. 
Datta  and  Jain  (1990)  [6]  investigated  the  response  of  an  articulated  tower  to  random  wave 
I  and  wind  forces.  In  the  derivation  of  the  single  degree  of  freedom  equation  of  motion  the  tower  is 
discretized  into  n  elements  having  appropriate  masses,  volumes  and  areas  lumped  at  the  nodes,  and 
there  is  viscous  damping.  The  equation  of  motion  is, 

I 

7(1  -b  ^(t))0  +  -h  72(1  +  =  F(t),  (7) 

where  is  the  time  varying  added  mass  term,  Jiu(i)  is  time  varying  buoyancy  moment  and  F(t) 

*  is  the  random  force  due  to  wave  and  wind.  The  Pierson-Moskowitz  spectrum  is  assumed  for  the 
wave  height  and  Davenport’s  spectrum  assumed  for  the  wind  velocity.  The  equation  is  solved  in 
the  frequency  domain  using  an  iterative  method,  which  requires  that  the  deflection  angle  6{t)  and 

I 

the  forcing  function  F{t)  be  decomposed  into  Fourier  series:  The  coefficient  of  the  sin  and  cos  are 
then  found  iteratively.  From  their  parametric  study,  they  concluded  the  following: 

'  1.  Nonlinearities  such  as  large  displacements  and  drag  force  do  not  influence  the  response  when 

only  wind  force  is  considered. 

2.  The  random  wind  forces  result  in  higher  responses  than  do  wave  forces. 

3.  The  r.m.s.  response  due  only  to  wind  forces  varies  in  a  linear  fashion  with  the  mean  wind 
velocity. 

In  a  later  paper,  Jain  and  Datta  (1991)  [8]  used  the  same  equation  and  the  same  method  of 
solution  to  investigate  the  response  due  to  random  wave  and  current  loading.  The  wave  loadings 
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(dr&g,  inertia  and  buoyancy)  are  evaluated  via  numerical  integration.  The  following  results  were 
obtained  from  the  parametric  study, 

1.  The  dynamic  response  is  very  small  since  its  fundamental  frequency  is  well  below  the  wave’s 
fundamental  frequency. 

2.  Nonlinear  effects  (drag  force,  variable  buoyancy)  have  considerable  influence  on  the  response. 

3.  Current  velocity  has  a  large  influence  on  the  response. 

Hanna  et  al.  (1983)  [15]  analyzed  the  nonlinear  dynamics  of  a  guyed  tower  platform.  The  tower 
is  represented  by  a  lumped  parameter  model  consisting  of  discrete  masses.  Each  mass  has  three 
degrees  of  freedom,  two  translations  and  one  rotation  about  the  vertical  axis.  The  external  forces 
on  the  structure  are  approximated  by  concentrated  forces  and  torques  at  the  nodal  points.  The 
equation  of  motion  is 

IMl(i)  +  (CKi)  +  (A-(u)|{«)  =  {i’ti.u.i)),  (8) 

where  [M]  is  the  total  mass  matrix  including  added  mass  terms,  [C]  is  the  structural  damping 
matrix  assumed  to  be  proportional  to  the  mass  matrix  and  [A^(u)]  is  the  total  nonlinear  stiffness 
matrix  that  includes  mooring  lines  effects,  soil  stiffness  and  geometric  stiffness.  {/*(<,«,«)}  is  the 
nonlinear  dynamic  load  vector  due  to  wave,  current  and  wind.  The  equation  is  then  solved  using 
normal  mode  superposition  and  the  response  is  calculated  at  each  time  step.  This  method  is  good 
only  if  the  nonlinearities  are  not  large.  Deterministic  and  random  loading  are  considered.  The 
solution  shows  insignificant  flexure  modes  while  the  torsional  one  has  a  noticeable  effect  on  the 
deck  rotational  response. 

Kanegaonkar  and  Haidar  (1988)  [10]  investigated  the  nonlinear  random  vibration  of  a  guyed 
tower.  They  included  nonlinearities  due  to  guylines  stiffnesses,  geometry,  load  and  material.  The 
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simplified  planner  equation  of  motion  is 


+  cd  +  =  M(0,  (9) 

where  if  is  a  spring  constant  depending  on  buoyancy,  gravity  and  guyline  horizontal  stiffness,  and 
ibi  is  a  constant  depending  on  the  guyline  vertical  stiffness.  M{t)  is  the  random  wave  loading.  The 
equation  is  then  solved  numerically  where  the  wave  height  is  defined  by  the  Pierson-Moskowitz 
spectrum.  It  was  seen  that  the  response  is  non-Gaussian  for  significant  wave  heights  greater  than 
5  m. 

Gerber  and  Engelbrecht  (1993)  investigated  the  response  of  an  articulated  mooring  tower  to 
irregular  seas.  It  is  an  extension  of  earlier  work  done  by  Thompson  et  al.  (1984)  [16].  The  tower  is 
modeled  as  a  bilinear  oscillator,  that  is,  a  linear  oscillator  with  different  stiffnesses  for  positive  and 
negative  deflections, 

mx  +  cx  +  {ki,k2)x  =  F{t).  (10) 

The  random  forcing  function  F{t)  is  assumed  to  be  the  sum  of  a  large  number  of  haurmonic  compo¬ 
nents,  each  in  different  frequency,  a  procedure  similar  to  that  proposed  by  Borgman  (1969)  [2].  The 
equation  is  then  solved  analytically  since  it  is  linear  for  each  half  cycle.  The  solution  is  obtained 
for  different  cases;  linear  oscillator  (both  stiffnesses  are  the  same),  bilinear  oscillator  and  for  the 
case  of  impact  oscillator  (a  rigid  cable)  in  which  oscillation  can  occur  only  in  one  half  of  the  cycle. 
For  future  study  they  suggest  to  include  nonlinecU'  stiffness  and/or  using  a  different  spectrum  to 
describe  the  wave  height. 
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1.3  Problem  Definition 


In  this  paper,  the  response  of  an  articulated  tower  submerged  in  the  ocean  is  investigated.  The  non¬ 
linear  differential  equation  of  motion  is  derived,  including  nonlinearities  due  to  geometry,  coulomb 
damping,  drag  force,  added  mass  buoyancy.  All  forces/moments  are  evaluated  analytically  ^nd 
explicitly  at  the  instantaneous  position  of  the  tower  and,  therefore,  they  are  time- depen  dent  and 
highly  nonlinear.  The  equation  is  then  numerically  solved  using  ‘ACSL’  -  Advanced  Continuous 
Simulation  Language  [1],  a  powerful  software  language  for  deterministic  emd  random  wave  loading 
using  the  Pierson-Moskowitz  wave  height  spectrum.  A  harmonic  and  subharmonic  solutions  for 
deterministic  wave  heights  are  obtained.  The  response  to  random  wave  heights  for  different  signif¬ 
icant  wave  heights  is  investigated,  the  influence  of  coulomb  damping  on  the  response  is  analyzed, 
and  chaotic  regimes  of  behavior  are  identified. 

The  distinctions  between  this  study  and  the  literature  with  which  we  are  aware  are  that, 

1.  A  sound  and  exact  derivation  of  the  nonlinear  equation  of  motion  is  provided. 

2.  All  terms  in  the  equation  of  motion  are  anal3rtically  derived. 

3.  Coulomb  friction  in  the  tower  hinge  is  added. 

4.  Usage  of  ‘ACSL’  for  the  numerical  solution  provides  an  easy  way  to  modify  parameters  and 
sensitivity  studies. 
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2  Problem  Description 


A  schematic  of  the  structure  is  shown  in  Fig.  1.  It  consists  of  a  tower  submerged  in  a  fluid  having  a 
concentrated  mass  at  the  top  and  one  degree  of  freedom  6  about  the  z  axis.  The  tower  is  subjected 
to  wave  loading.  Two  coordinate  systems  are  used;  one  fixed  (x,i/,z)  and  the  second  attached  to  the 
tower  (x',y',z').  All  forces/moment  are  derived  in  the  fixed  coordinate  system,  which  means  that 
the  tower  rectilinear  velocity  is  resolved  into  x,  y  coordinates.  The  motion  of  the  tower  is  assumed 
to  be  only  in  plane. 

This  problem  has  similarities  to  that  of  an  inverted  pendulum,  but  due  to  the  presence  of  gravity 
waves,  additional  considerations  are  made  : 

1.  Forces  due  to  buoyancy  and  vertical  wave  velocity  are  summed  and  denoted  as  To. 

2.  Drag  forces  proportional  to  the  square  of  the  relative  velocity  between  the  fluid  and  the  tower 
need  to  be  considered. 

3.  Fluid  inertia  forces  due  to  fluid  acceleration  are  part  of  the  loading  environment. 

4.  Fluid  added  mass  is  directly  included  in  the  inertia  forces. 
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3  Equation  of  Motion 


The  equation  of  motion  is  derived  using  Newton’s  second  law,  setting  all  applied  moments  acting 
on  the  tower  about  hinge  equal  to  the  dynamic  moment.  We  will  find  the  equation  of  motion  to  be 
of  the  form, 

J|j(«)l»  =  Ev  [(.w, /(»),%»(*)] ,  (11) 

where,  J\g{9)]  is  the  effective  position -dependent  moment  of  inertia,  g{0)  and  f{6)  are  nonlinear 
functions  of  6  (trigonometric  functions),  u>  is  the  wave  frequency,  and  YIM  is  the  sum  of  all  external 
moments  that  act  on  the  tower.  Certain  assumptions  have  been  made  in  deriving  the  nonlinear 
equation  of  motion.  These  are  listed  below. 


3.1  Assumptions 

1.  The  tower  stiffness  is  infinite;  El  =  oo. 

2.  The  hinge  has  coulomb  damping. 

3.  The  tower  has  a  uniform  mass  per  unit  length,  m,  and  is  of  length  /  and  diameter  D. 

4.  The  tower  is  a  smooth  slender  structure  with  uniform  cross  section. 

5.  The  end  mass  M  is  considered  to  be  concentrated  at  the  end  of  the  tower.  (It  has  no  volume.) 

6.  The  tower  length  is  greater  than  the  fluid  depth,  but  the  dynamics  is  not  limited  to  the  case 
of  the  mass  always  being  above  the  mean  water  level. 

7.  The  structure  is  statically  stable  due  to  the  buoyancy  force. 

8.  The  waves  are  linear  having  random  height. 

9.  Morison’s  fluid  force  coeflScients  Cd  and  Cm  are  constant. 

10.  The  center  of  mass  (c.g.)  of  the  tower  is  at  its  geometric  center. 

11.  Currents,  wind  and  wave  slanuning  forces  are  not  included. 

3.2  Forces/Moments  acting  on  the  tower 

Fig.  2  describes  the  external  forces  acting  on  the  tower.  These  are  : 

1.  To  is  a  vertical  buoyancy  force. 

2.  Fy,  Fh  are  the  vertical  and  horizontal  fluid  forces  due  to  fluid  drag  and  inertia. 

3.  Mg,  mlg  are  the  forces  due  gravity. 

We  next  describe  these  forces  and  moments. 


X 


Figure  2:  Ebctemal  Forces  on  the  Tower 


3.2.1  Inertia  Moment 


The  inertia  moment  equals  the  total  moment  of  inertia  of  the  tower  multiplied  by  the  angular 
acceleration  of  the  tower  plus  the  fluid  added  mass  term, 


(12) 


where  Jo  is  the  moment  of  inertia  of  the  tower  about  point  ‘o’  and  Mji  is  the  fluid  added  mass 
moment.  The  tower  moment  of  inertia  is 

Jo  =  im;=.  (13) 

The  added  mass  force  per  unit  length  is 

F,i  =  (14) 

where  Ca  is  the  added  mass  coefficient  which  equals  Ca  =  Cm  —  1,  Cm  is  the  inertia  coefficient 
and  V  is  the  tower  acceleration.  In  order  to  find  the  added  mass  term  the  tower’s  acceleration  is 
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derived  by  taking  the  second  time-derivative  of  the  tower  radius  vector  in  the  x,  y  coordinate, 
R  =  x'  cos  9x  -f-  x'  sin  By 

JTJ 

-7-  =  V  =  —x'O  sin  Ox  -f  x'B  cos  By 
at 

=  V  =  — x'(dsind -I- 0^cos0)x -f  x'(dcos0  —  (j^sin0)y.  ( 

at 


Replacing  x'  =  ■—  we  find 


R  =  XX  -1-  X  tan  By 

V  =  —xB  tan  Bx  -f  xBy 

V  =  — x(d  tan  0 -I- 5*)x -i- x(^  —  tan  )y 


and  the  fluid  added  mass  moment  is 


Mfi  —  j  R  X  Ffidx. 


Substituting  the  expression  for  V  in  Fji  and  evaluating  the  cross  product  yields, 


x^(l  -t-  tan*^)dx 


where  L  is  the  projection,  in  the  x  direction,  of  the  submerged  part  of  the  tower,  and  is  quantified 
later.  Thus,  the  fluid  added  mass  moment  is 

Mji  =  \CaP1^—L\\  +  (19) 

and,  assuming  that  the  fluid  added  mass  due  to  the  end  mass  is  negligible,  the  total  inertia  moment 


Mi  =  Je//0, 


where  Je//»  the  total  moment  of  inertia,  is 


Jfi  =  \{ml  -1-  M)fi  -i-  1,^(1  -H  tan^B). 

3  3  4 
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3.2.2  Moments  due  to  Gravity  and  Buoyancy 

The  following  moment  is  due  to  terms  that  do  not  depend  on  the  fluid  velocity,  such  as  gravity  and 


buoyancy  forces, 


Mg  =  Toll,  —  {M  +  rn^)gl  sin  0, 


where  g  is  the  gravitational  constant.  To  is  the  time  dependent  buoyancy  force  and  4  is  its 


To  =  pgVo 

where  Vo  is  the  volume  of  the  submerged  part  of  the  tower.  Assuming  that  Z>  «  d  we  find 

To  =  pgif—L,, 

4 

and  L„  which  is  the  length  of  the  submerged  part  of  the  tower,  is 

r  d  + 

’  cos9 


d  is  the  mean  water  level,  p  is  the  fluid  density,  is  the  wave  height  elevation  given  by 


i)  =  2^  cos{ky  -  wt  +  e), 


where 


=  gkiaxih{kd) 


H  is  the  wave  height,  u  is  the  wave  frequency,  and  k  is  the  wave  number.  Since  we  are  interested 
in  7  at  the  instantaneous  position  of  the  tower  and  at  x  =  d  with  y  =  d  tan  9, 


T}{9,t)  —  -H  cos{kdiaLn9  —  u;t  +  e). 


15 


The  buoyant  force  acts  at  the  center  of  mass  of  the  submerged  part  of  the  tower.  If  we  consider  the 
tower  to  be  of  cylindrical  cross-section  then  the  center  of  mass  in  x\  y'  coordinates  is 


= 


It 


1 ,  2. 


(29) 


Transforming  to  x,  y  coordinates  we  '^.nd 


4  =  /j  cos  6  +  if  sin  6 

1 

4  =  tan^  ^  cos  ^ -h -I,, -i- —  tan^^sing. 

x6.uj  2 


(30) 


3.2.3  Moriflon’s  Equation  for  Wave  Force 


In  general,  the  fluid  forces  acting  on  a  slender  tower  are  of  two  types:  drag  and  inertia.  The  drag 
force  is  proportional  to  the  square  of  the  relative  velocity  between  the  fluid  and  the  tower,  and  the 
inertia  force  is  proportional  to  the  fluid  acceleration.  The  drag  and  inertia  forces  per  unit  length 
are  approximated  by  Morison’s  equation, 

F/,  =  CdP^  1  U  -  V  I  (U  -  V)  -1-  (31) 


where  F//  is  the  fluid  force  per  unit  length,  and  U,  V  are  the  fluid  and  tower  absolute  velocities, 
respectively.  In  order  to  find  the  moments  due  to  the  fluid  forces,  the  fluid  and  tower  absolute 
velocities  are  divided  into  two  components:  u  in  the  x  direction,  w  in  the  y  direction.  Then  the 
horizontal  and  vertical  forces  are  evaluated. 

Assuming  a  linear  deterministic  wave  field  initially,  the  horizontal  and  vertical  fluid  velocities 
are  respectively  (Wilson  [17]  pp.  84)  : 


u 


1  cosh  kx 

2  ^  sinh  kd 


cos{ky  —ut) 
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w 


(32) 


1  „  sinh^x 

2  sinh  kd 


sm{ky  —  u;<), 


and  the  respective  accelerations: 


u 

w 


1„  2sinhA;x  ,,  , 

-2"“  riSw 


(33) 


To  evaluate  the  fluid  forces,  the  velocities  and  accelerations  have  to  be  defined  along  the  tower. 
Substituting  this  expression  for  y  =  x  tan  6  in  equation  (32)  yields  the  fluid  horizontal  and  vertical 
velocities  along  the  tower, 


u 

tv 


\h. 


cosh  A;x 
sinh  kd 
sinh  kx 
sinh  kd 


cos(ibx  tand  —  ujt) 
sin(fcxtan0  —  ut). 


(34) 


Fig.  3  depicts  sample  velocity  profiles  of  the  fluid  and  the  tower  (equation  (16))  along  the  tower 
length.  The  velocities  are  time  dependent  and  in  the  figure  they  are  shown  at  a  particular  instant 


of  time. 


Figure  3:  Tower  end  Fluid  Velocity  Profiles  (Solid  Line  =  Fluid,  Dashed  Line=  Tower) 


Takiog  the  derivative  of  the  velocities  with  respect  to  time,  noticing  that  0  is  time  dependent, 
leads  to  the  accelerations 

«  =  ( —u  -h  6  sin(A;x  tan  6  — 

2  \  cos*  V )  smh  kd 

1  rr  /  n  \sinhA:x  ,,  . 

w  =  -/fw  1  — w  +  — rr  1  cos(fcj tan ^  —  u>t).  (35) 

2  V  cos*  smh  kd  ^  ^  ’ 

The  moments  due  to  the  fluid  velocity  and  acceleration  are  as  follows.  The  inertia  moments  are 

A/m  =  Cmpk--  /  uxdx 
4  Jo 

Mtv  =  Cm  fur— ^  /  wxtan^dx,  (36) 

4  JQ 

and  the  drag  moments  are  given  by 

D 

Moh  -  ^DP-^  I  ti  ~  t;„  I  (u  -  Vy)xdx 

D 

Mdv  =  CdP~t  I  I  u>  ~  V*  I  (tn  —  Vx)xtanxOdx.  (37) 

2  Jo 

where  Mmi  A/m  are  the  inertia  moments  and  A/pfc,  Mdv  are  the  drag  moments  due  to  horizontal 
and  vertical  accelerations,  respectively.  The  moment  arm  for  the  vertical  force  is  y  =  x  tan  6.  L, 
which  is  the  upper  limit  of  the  integral,  depends  on  the  angle  9  as  follows  ; 

I  cos  9  if  d  >  /  cos  9 

L  =  (38) 

d  if  d  <  /  cos  9. 

This  means  that  the  moments  are  calculated  only  for  the  submerged  part  of  the  tower. 

Substituting  the  velocities  (34)  and  accelerations  (35)  into  the  moment  equations  yields  the 
inertia  moments 


A/m  =  CmP^^  f  \hu)  f-cj  +  9  — ;  cos(A:j  tan  g  -  ut)  xtan^dx,  (39) 

4  do  [2  \  co8^9J  sinhkd 
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and  the  drag  moments 


Mok  =  CdP^  f  ( f  j  cos{kx  tan 6  —  u>t)  --  x6\  xdx 
2  Jo  \2  sinhica  / 

^Dv  =  CdP^  f  i sm(kx  tan  6  —  ut)  +  x6  tan  ^  )  xt&n  Odx. 
2  Jo  \2  smhKa  J 


The  direction  of  the  drag  moment  depends  on  the  sign  of  the  relative  velocity  between  the  fluid 
and  the  tower.  In  equation  (40),  it  is  assumed  that  this  velocity  is  positive.  When  solving  the 
governing  nonlinear  differential  equation  numerically,  the  relative  velocity  must  be  evaluated  for 
each  time  step  to  determine  that  sign. 

The  integrals  in  equations  (39)  and  (40)  are  very  complicated  to  evaluate  analytically.  Although 
this  work  is  concerned  with  large  angles,  we  are  interested  in  examining  the  limitations  of  the  ‘small 
angle’  assumption.  To  do  this,  analytic  expressions  for  the  above  integrals  are  evaluated  using 
‘MAPLE’  about  y  =  0,  and  compared  with  numerical  integrations.  This  exercise  shows  that  this 
assumption  is  very  good  and  the  error  for  0  <  30®  is  less  the  1  %.  Thus,  we  have  at  our  disposal 
accurate  analytical  approximations. 

Substituting  (34),  (35)  into  the  moment  formulas  and  integrating  leads  to  the  inertia  and  drag 
moments,  respectively, 


=  +  kL  +  e"'"  +  1  -  26“-) 


(2  -  2  -  2  -2kL  +  0 


cos^d 


- 1)  “ 

(2  4-2-2  +  2  jfcl  -1-  -  4  0 1 

cos^  0  J  ’ 
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_  Die^L*  /(-jb  +  2Pl)c*^^  2ib  +  4PL2  -2kn-k\^^ 
Mo.  =  CoPjl^—^\^ - +  + 


'  (32  fci  -  32  -  16  k^L^)  c*^  32  ibi  +  32  +  16  k^L^ 

16  P  16 


)«Bi|, 


Mov  —  Cop 


2\  4  ( 


(-k  +  2k^L)e^‘-^  2k-Ak^P  -2k^L-k 


- ^  A.X/  -  '''1^0  >* 

+  16t3e‘‘=  r“ 


'(32itI-32-16PI.“)e“  4  -li k^ -  32  kL- 


+  F  + 
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i?)« 


tan*  5^2, 


where  Bi,  are  defined  as 


B,  =  . 

"  2  sinMJbd) 

The  drag  moment’s  direction  is  opposite  to  the  relative  velocity  between  the  fluid  and  the  tower. 
Since  this  velocity  depends  on  depth  (z  axis)  and  may  have  different  signs  along  the  length,  we 
will  find  its  average,  and  the  sign  of  the  average  relative  velocity  will  determine  the  direction  of  the 
resultant  drag  moment. 

The  general  expression  for  the  average  relative  velocity  is  given  by  : 

—  1 


Vre/  =  J  [  (U  -  y)dx. 
mj  Jo 


We  have  to  evaluate  the  horizontal  V^.,  and  vertical  relative  velocities.  Thus,  we  have 


=-r(- 

lJo  \2 


cosh  lex  •  \ 

Hu)  .  ,  ,  ,  cos{ibz  tan  9  —  ut)  —  x9 )  di, 
smh  kd  ) 


and  after  integrating  we  get 


— u  9L  \Hu}co^^9  ,  ,  .  ,,  ,  «  ,, 

.in'h(M)  (”°h(tl)c°s(*l.t»n(l  -wt) 

+  cosh(fci)tan0sin(A:Z<tan0  —  u><)  +  2  tan^  sin(ci;t)}  . 
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Similarly 


LU  V2 


„  smhkx  .  , 

nu)  .  .  .  .  sm( Kx  tan  0  —  ut)  +  x0  tan  0  dx, 
Sinn  kd  ) 


and  after  integrating 


^  0Lia.xi0  \H(jJcos^0  .  ,,  ,  .  , 

^rei  —  - - \- -f-r-r-rjj-r  {—smh{kL)taxi0  cjos(kLtaa0  —  ut) 

2*  wCAj  SlUu^flbUy 

+  cosh{kL)  sin(kLtaii0  —  u/t)  +  2  sin(u;t)}  . 


When  we  subsequently  solve  the  equation  of  motion  numerically,  we  will  check  the  sign  of  the 
relative  velocities  to  set  the  correct  direction  of  the  drag  moments  and  also  to  check  the  sign  of 
(d  —  I  cos  0)  to  set  the  limits  of  integration. 


3.2.4  Damping  Moment 

The  tower  hinge  is  assumed  to  dissipate  energy  via  coulomb  friction.  In  this  section,  this  fric¬ 
tion/damping  moment  is  evaluated.  The  damping  force  is  equal  to  the  product  of  the  normal  force 
N  and  the  coefficient  of  friction  /i.  It  is  assumed  to  be  independent  of  the  velocity,  once  the  motion 
is  initiated.  Since  the  sign  of  the  damping  force  is  always  opposite  to  that  of  the  velocity,  the 
differential  equation  of  motion  for  each  sign  is  valid  only  for  a  half  cycle  interval.  The  friction  force 
is 

Ffr  =  Nftlsgn(0)].  (47) 

The  normal  force  is 

N  =  ^FgCos0 +  ^FySm0,  (48) 

where  ^Fx,  Fy  are  the  total  forces  in  the  x,  y  directions,  respectively.  These  forces  are  due  to 
gravity,  buoyancy  and  fluid  drag  and  inertia, 

Y^Fx  =  ro-/i  +  F/„  +  Fx,„ 
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=  ^//i  +  ^DA, 


(49) 


where  To  is  the  buoyancy  force  given  in  equation  (24),  Fg  is  the  gravitational  force, 


Fg  =  (ml-i-M)g,  (50) 

and  Fohi  T'dv,  Flh,  Fiv  are  the  fluid  drag  and  inertia  forces  in  the  horizontal  (y)  and  vertical  (x) 
directions,  respectively.  These  forces  are  calculated  using  Morison’s  equation  in  a  similar  way  as 
the  moments  and  are  given  by, 

Fdh  -  <^d/>2  I  3  8ifce“7  J 


Fjh  =  Cmpv- 


{((  e'‘^^kL  -kL-  - 


{(( 


2 


- + k-^ 


08^0  2  c*^Ar 


)'■) 


Fiv  -  2ibV^  cos^O  2c*H  ‘^Jb] 

If  we  assume  a  hinge  radius  Rt^,  then  the  damping  moment  is 

Mjr  =  i?A  (H  cos  ^  F,  sin  ff)  fil^gn(0)J.  (52) 

Finally,  adding  all  the  external  moments,  equations  (22),  (41),  (42),  (52),  and  setting  them  equal 
to  the  total  inertia  moment,  equation  (20),  yields  the  governing  nonlinear  differential  equation  of 


motion. 


Jef/0  =  —M/r  —  Mg  +  Af/fc  —  Mjy  +  Moh  ~  Mdv' 


Both  sides  of  this  equation  are  nonlinear  functions  of  0,  0,  t  and  u;  : 

Jejf(0)6  =  -Mfr(0,u,t,sgn(6))- Mg{0,u;,t)  + - 

Miv{0,O,(^,t)  +  MDk{0,i*f,t)  -  MDv(0,9,iJ,t), 
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which  agrees  with  the  general  form  of  equation  (11). 


4  Numerical  Solution 

The  numerical  solution  for  the  nonlinear  differential  equation  (53)  was  performed  using  ‘ACSL’ 
and  the  analyses  of  the  results  was  performed  using  ‘MATLAB’.  The  ‘ACSL’  code  written  for  this 
application  is  available  upon  request.  In  order  to  build  confidence  in  the  solution,  some  test  cases 
have  been  solved  using  the  following  physical  parameters  : 

Tower  properties 

1.  /  -  Length  of  the  tower  =  400  (m) 

2.  D  -  Tower  diameter  =  15  (m) 

Z.  M  -  End  mass  =  2.5  *  lO’  (Kg) 

4.  m  -  Tower  mass  per  unit  length  =  2  ♦  10®  (^) 

5.  /i  -  Friction  coefficient  =  0. 1-0.4 

6.  Rh  -  Hinge  radius  =  1.5  (m) 

Fluid  properties 

d  -  Mean  water  level  =  350  (m) 

2.  Cd  -  Drag  coefficient  =0.6  -  1.0 


3.  Cm  -  Inertia  CoeflBcient  =1.5 

4.  p  -  Water  density  =  1025  (^) 

5.  w  -  Wave  frequency  =  0.2  -  1.0  (^) 
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4.1  Response  for  Deterministic  Wave  Heights 


In  this  section  the  response  of  the  tower  to  deterministic  wave  heights  is  established.  Free  vibration 
and  stability  analyses  are  performed. 

First,  the  natural  frequency  of  the  tower  is  found  as  if  it  was  an  upright  pendulum  subjected 
to  a  constant  tension  force  To  and  gravitational  force  Mg.  Fig.  4  shows  the  tower  response  to  an 
impulse  for  Cd  =  0  and  H  =  0  \n  the  time  and  frequency  domain.  The  response  is  harmonic  in 
the  tower  natural  frequency  =  0.026  (Hz),  and  it  agrees  with  the  calculation  for  the  natural 


frequency  of  a  pendulum  : 


fin  = 


i.  Tp/fc  -  (0.5m/  -b  M)gl 
2ir  \  J,ii 


(55) 


Tkiit  Donah  RMponn  ID 


ThalOM) 


Ho.'uonoyl^'l)) 


Figure  4:  Tower  Natural  Frequency,  Co  =  0 
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Fig  5  (a,b)  shows  the  same  but  for  /i  =  0  and  C©  =  1,  the  decay  here  is  not  linear  since  the 
drag  force  is  propotional  to  the  velocity  squared.  Fig  5  (c,d)  shows  the  free  vibration  with  frictional 
damping,  /i  =  0.1,  in  the  time  and  frequency  domain.  Here  the  amplitude  decays  linearly  with 
time  as  expected  when  coulomb  damping  is  present.  From  the  figure  we  can  evaluate  the  equivalent 
damping  ratio  for  Cd  =  0  to  be  ^  =  0.02.  Since  the  damping  is  nonlinear,  the  ‘natural  frequency’ 
smd  its  multipliers  are  seen  in  the  figure. 

TneDocnain  Response  (a)  Tim  Domain  Response  (c) 
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Figure  5:  Tower  FVee  Vibration  with  Damping,  (a,  b)  ^  =  0.1,  (c,  d)  Cxj  =  1 


Fig.  6  shows  the  response  in  the  time  domain  and  frequency  domain  for  /f  =  1  (m)  and  wave 
loading  frequency  w  =  0.064  {Hz).  Figures  (a,  b)  are  with  /i  =  0  and  (c,  d)  are  with  n  =  0.4.  In  the 


frequency  domain,  the  tower  ‘natural  frequency’  and  the  wave  natural  frequency  are  clearly  seen. 
It  can  be  seen  that  the  friction  has  a  damping  and  stabilizing  effect  on  the  system. 


Figure  6:  Tower  response  to  wave  excitation,  Time  and  Frequency  Domain,  H  =  l  (m) 


Since  this  is  a  nonlinear  system,  its  response  to  harmonic  excitation  at  the  system’s  'natural 
frequency’  or  its  multipliers  can  be  multivalued,  indicating  the  occurrence  of  a  ‘jump’.  Fig.  7  shows 
the  tower  response  to  harmonic  wave  excitation  at  u;  =  n„  (harmonic  solution)  with  and  without 
damping.  If  the  system  was  linear,  such  loading  would  have  caused  resonance  instability.  But  since 
the  system  is  nonlinear,  an  amplitude  change  (beating)  can  be  seen  Fig.  7  (a),  indicating  that 
‘jumps’  occur  from  one  amplitude  to  the  other  (Wilson  [17]  pp.  147).  It  can  also  be  seen  that 
although  the  excitation  frequency  is  constant,  the  frequency  response  is  not  Fig.  7  (b).  Again  the 
reason  is  the  nonlinear  system  characteristics.  From  the  response  with  friction  (figures  (c,  d))  we 
see  that  the  response  is  not  smaller  since  it  is  unstable,  but  the  ratio  between  the  large  and  small 
amplitude  of  the  beating  gets  smaller. 

TmOofflMRHponMM  TiMDonuiRMpamjc) 


Oil 


- - - 1 

0  500  1000  1500 

TimM 


0.4r 


•0.4' - -  - - -I 

0  500  1000  1500 

Tim(mc) 


Figure  7:  Tower  response  to  harmonic  wave  excitation  at  the  ‘Natural  Frequency’  -  Beating  Phenomenon 


Fig.  8  shows  the  response  due  to  haxmonic  wave  excitation  at  twice  the  system  ‘natural  fre¬ 
quency’  (sub-harmonic  solution),  with  (c,  d)  and  without  friction  (a,  b).  We  can  see  the  beating 
phenomenon,  and  from  the  frequency  domain  we  see  that  the  response  consists  of  the  ‘natural 
frequency’and  it  multipliers.  This  is  due  to  the  nonlinear  nature  of  the  system. 
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Figure  8:  Tower  response  to  sub-harmonic  wave  excitation  at  twice  ‘Natural  Frequency’ 
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4.1.1  Chaotic  Response 


Because  of  the  nonlinear  characteristics  of  the  system,  the  governing  differential  equation  of  motion, 
chaotic  motion  can  occiire  under  certain  conditions.  Fig.  10  shows  the  response  of  the  tower  to  a 
harmonic  wave  excitation  at  a  frequency  u  =  0.95n„.  We  see  that  the  response  ‘jumps’  between 
two  stable  equilibrium  states.  This  phenomenon  can  imply  a  chaotic  system. 
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Figure  10:  Tower  response  to  harmonic  wave  excitation  at  the  ‘Natural  Frequency’,  Cd  =  0.6 
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Moon  [11]  provides  several  ways  to  identify  chaotic  vibrations.  We  next  use  three  ways  to  prove 
that  the  system  under  consideration  is  chaotic.  Fig.  11  describes  the  response  of  the  tower  to  a 
harmonic  wave  excitation  having  a  frequency  of  a;  =  0.09  {Hz).  Fig  11  (a,  b)  shows  the  response 
in  the  time  and  frequency  domains  and  the  phase-plane  trajectory  (c).  From  the  time  history  it  is 
clearly  seen  that  the  system  jumps  between  two  stable  equilibrium  states.  The  frequency  domain 
shows  multiharmonic  energy  although  the  system  has  one  degree  of  freedom;  another  indicator  of 
chaos.  Finally,  from  the  phase-plane  trajectory  we  see  that  the  orbits  never  close  or  repeat.  Thus, 
we  conclude  that  this  model  has  the  characteristics  of  a  chaotic  system.  This  chaotic  phenomenon 
appeared  only  when  the  wave  height  was  larger  than  10  (m). 
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Figure  11:  Tower  response  to  a  harmonic  wave  excitation,  u  =  0.09  {Hz)  -  Chaos 


4.2  Random  Wave  Height 


In  this  section  the  tower  response  to  random  wave  height  excitation  is  investigated.  The  wave 
height  distribution  is  generally  expressed  in  the  form  of  a  power  spectral  density.  For  a  simulation 
of  the  response  in  the  time  domain,  the  power  spectra  of  the  wave  height  is  transformed  into  a  time 
history.  This  is  accomplished  using  a  method  by  Borgman  [2],  described  in  Wibon’s  book  [17],  and 
presented  here. 


4.2.1  Synthesis  of  Time  Histories  firom  Power  Spectra 


The  wave  elevation  Tf{x,t)  can  be  expressed  as 


i/(x,f)  =  J°°  co8{kx  —  wt  +  e)y/ A^{u>)du), 


(56) 


where  A‘^{oj)  is  the  amplitude  spectnun  (height)  and  e  is  a  random  phase  angle  having  a  uniform 
distribution  over  an  interval  0  to  2ir.  To  evaluate  the  integral,  the  spectrum  is  discretized  into  equal 
areas  (not  equal  frequencies).  This  procedure  avoids  the  presence  of  periodicities  in  the  resulting 
time  history.  Consider  the  following  partition  ; 


Ck\)  <  <  W2  ^  ••••  ^  ~ 


(57) 


where  uo  is  a,  small  positive  value  and  F  is  the  frequency  above  which  the  spectral  amplitude  is 
practically  zero.  Let 


Au;„  =  u;„  —  u)n-i 

Un  -  (jJn-1 


U>„  = 


;  n=l,2, 


(58) 

(59) 
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where  N  is  the  number  of  partitions  of  the  spectra.  Now  the  integral  can  be  approximated  as  a 


vix,  0  =  iz  COs(i„X  -  £n)y/A^{Un)^U}n,  (60) 

nsl 

where  a>*  =  k^g  and  ()  represents  the  average  of  each  parameter.  Let  5(u;„)  represent  the  cumulative 
area  under  the  spectral  density  curve,  or 


Thus, 


5(w„)  =  5Z  ^^(w„)Au?„. 


«  5(u;„)  -  5(t4>n-i)  = 


where  is  a  constant  to  be  determined.  It  follows  that 


Na^  =  «  <S(oo)  s=  /  A^{ut)duf. 

Jo 


The  Pierson-Moskowitz  spectrum  for  the  wave  height  is 


where  Ao  and  B  are  constants  defined  by 


Ao  =  8.1  X  10"V 


B  = 


Then,  from  equations  (61)  and  (64), 


and,  therefore. 


4B 

4BN' 
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Now  the  partition  frequencies  can  be  determined  for  w  =  F, 


Because  the  partition  is  of  equal  areas 


5K)  =  ^S(F)  =  =  5(F)«»'®*e-®'“‘-, 


it  foUows  that 


£L^B/F*  ^  gB/wJ 


Solving  equation  (71)  for  u;„  leads  to  the  partition  frequencies 


-  (  ^  V'”  n  -  1  2  N 

~\\n{N/n)  +  B/F*)  ’ 


Therefore,  the  wave  elevation  ri{x,  <)  can  be  approximated  by 

N 

fj{x,  <)  =  <1 53  +  £«)•  (73) 

nsl 

The  wave  loading  on  the  tower  is  a  function  of  wave  velocity  and  acceleration.  Therefore,  in 
our  numerical  studies,  these  have  to  be  expressed  as  fimctions  of  the  approximate  wave  elevation. 
Thus,  in  the  expressions  for  wave  vdodty  and  acceleration,  the  following  substitutions  are  made 

H  =>  (74) 

SABN  ^  ’ 


where  n  =  l,2,...iV.  For  example,  the  horizontal  velocity  and  acceleration  will  become 

>4o  cosh /l  *  a 

”  =  5  4Bjv“’riSW 


i  =  E- 


ABN 
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As  mentioned  earlier,  the  wave  height  is  assumed  to  have  the  Pierson- Moskowitz  spectrum 


Sr.io,)  = 


0.0081  * 
- 5 - 

/.l® 


(76) 


The  significant  wave  height,  H,,  changes  the  maximum  wave  height  and  the  wave  frequency  distri¬ 
bution  as  shown  in  Fig.  12  for  H,  =  15  and  9  (m). 


Figure  12:  Examples  of  Two  Pierson-Moskowitz  Wave  Height  Spectra 
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4.2.2  Reaponse  for  Random  Wave  Height 


In  this  section,  the  influance  of  different  significant  wave  heights,  drag  coefficient  ,  hinge  friction 
and  non-zero  intial  condition  on  the  response  is  investigated.  Fig.  13  compares  the  tower  response 
for  /f,  =  9  and  15  (m).  Here  we  enlarged  the  buoyancy  force  so  that  the  tower  ‘natural  frequency’ 
is  Hn  =  0.07  {Hz)  ~  0.44  {radjsec).  It  can  be  seen  that  the  deflection  angle  for  H,  =  9  (m)  is  of 
the  same  magnitude  as  for  for  H,  =  15  (m),  although  the  later  maximum  wave  height  is  four  times 
larger  then  the  former  (as  can  be  seen  from  Fig.  12).  The  reason  is  that  the  natural  frequency  of 
the  tower  coinsides  with  where  most  of  the  spectral  aiergy  for  H,  =  9  (m)  is  located,  as  shown  in 
Fig.  12.  Since  the  lowest  wave  freqency  is  about  u  =  0.2  (radlsec)  articulated  tower  are  designed 
to  have  a  ‘natural  frequency’  lower  than  that. 
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Figure  13:  Tower  Response  -  Deflectitm  angle  and  Total  Moment  for  H,  =  15  (dashed  line)  and  9  (solid  line)  (m) 
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A  comparison  of  the  tower  response  for  different  values  of  the  drag  coefficient  are  plotted  in 
Fig.  14.  This  figure  shows  the  deflection  angle  and  the  total  drag  moment,  A/j  =  Moh  -  Mdv,  for 
Cd  =  0.6  and  1,  and  H,  =  9  (m).  It  can  be  seen  from  the  figure  that  the  deflection  angle  and  the 
drag  moment  are  larger  for  Cd  =  1>  since  the  drag  moment  in  the  Morison  equation  is  proportional 
to  the  drag  coefficient.  The  results  were  similar  for  different  significant  wave  heights. 


Tine  (sec) 


Figure  14:  Tower  Response  -  Deflection  angle  and  Total  Drag  Moment  for  H,  =9  (m)  and  Cd  =  0.6  (dashed  line) 
and  1  (  solid  line). 
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The  friction  effect  on  the  response  is  shown  cleaxly  in  Fig.  15.  Here  the  significant  wave  height 
is  H,  =  9  and  Cn  =  1.  It  is  seen  that  the  response  is  much  smaller  and  smoother  due  to  the 
additional  friction  in  the  system. 
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Figure  15:  Tower  Response  -  Deflection  angle  -  Time  and  Frequency  Domain  for  H,  =  9  (m)  and  Cc  •=  1 ,  (a,  b) 
/I  =  0  and  (c,  d)  /<  =  0.4 
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The  response  for  nonzero  initial  conditions  is  depicted  in  Fig.  16.  The  response  for  =  0)  = 
0.01(rad/scc)  (a,  b),  and  6{t  =  0)  =  0.05(rad)  (c,  d).  The  drag  coeiBcient  is  Cd  =  0.6  and  H,  =  9 
(m). 
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Figure  16:  Tower  Response  -  (a,  b)  Deflection  angle  9{t  =  0)  =  0.05  (rad)  and  (c,  d)  0{t  =  0)  —  0.01  {rad/sec). 


5  Discussion  and  Summary 


The  nonlineau*  differential  equation  of  motion  for  an  articulated  tower  submerged  in  the  ocean  is 
derived.  Geometric  as  well  as  force  nonlinearities  are  included  in  the  derivation.  The  wave  velocities 
and  accelerations  are  determined  at  the  instantaneous  position  of  the  tower,  a  fact  that  added  to  the 
nonlinearities  of  the  equation.  The  equation  is  solved  numerically  using  ‘ACSL’  for  deterministic 
and  random  wave  loading. 

The  response  of  the  tower  to  harmonic  wave  excitation  at  its  ‘natural  frequency’  and  at  twice  its 
‘natural  frequency’  demonstrates  beating,  where  the  amplitude  varies  between  two  extremes.  This 
beating  is  due  to  the  nonlinear  behaviour  of  the  system.  Coulomb  damping  reduces  the  beating 
phenomenon  and  the  response  amplitude,  so  it  has  a  stabilizing  effect  on  the  system.  When  the 
system  is  excited  at  an  arbitrary  frequency,  and  the  wave  height  is  greater  than  about  10  (m),  the 
response  ‘jumps’  between  two  stable  equilibria,  exhibiting  chaotic  behaviour. 

To  solve  the  equation  for  random  wave  loading,  the  Pierson-Moskowitz  spectrum  that  describes 
the  wave  height  distribution  was  first  transformed  into  a  time  history.  The  equation  was  solved  for 
two  significant  wave  heights.  Again  the  response  was  periodic  consisting  of  the  fundamental  tower 
frequency  and  its  multipliers.  For  significant  wave  h«ght  of  9  (m),  the  response  was  larger  than  that 
for  15  (m),  since  in  the  former  the  tower  ‘nattiral  frequency’  coincides  with  the  frequencies  where 
most  of  the  energy  is  located.  The  response  with  coulomb  damping  shows  that  friction  stabilizes 
the  system.  Notice  that  in  order  to  reduce  stresses  in  the  structure,  the  friction  moment  has  to  be 
low  enough  so  that  the  tower  can  comply  with  the  wave  loading. 

A  more  realistic  model  having  two  angular  degrees  of  freedom  is  being  analyzed  at  the  present 
time.  The  response  due  to  wave,  current  (colinear  and  non-colinear)  and  vortex  shedding  loading 
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is  investigated  loading  and  results  will  be  published  in  the  near  future. 
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ABSTRACT 


This  paper  presents  studies  on  the  response  of  an  ar¬ 
ticulated  tower  in  the  ocean  subjected  to  deterministic 
and  random  wave  loading.  The  tower  is  modeled  as 
an  upright  rigid  pendulum  with  a  concentrated  mass 
at  the  top  and  having  one  angular  degree  ci  freedom 
about  a  Unge  with  coulomb  damping.  In  the  derivar 
tion  of  the  differential  equation  of  motion,  nonlinear 
terms  due  to  geometric  (large  angle)  and  fluid  forces 
(drag  and  inertia)  are  included.  The  wave  loading  is 
derived  using  Morison’s  equation  in  which  the  velocity 
and  acceleration  of  the  fluid  are  determined  along  the 
instantaneous  position  of  the  tower,  causing  the  equar 
tion  of  motion  to  be  highly  nmilinear.  F^irthermore, 
since  the  differential  equation’s  coefficients  are  time- 
dependent  (periodic),  parametric  instability  can  occur 
depending  on  system  parameters  such  as  wave  height 
and  frequency,  buoyancy,  and  drag  coeflicient.  The 
nrailinear  differential  equation  is  then  solved  numeri¬ 
cally  using  ‘ACSL’  software.  The  response  of  the  tower 
to  deterministic  wave  loading  is  investigated  and  a  star 
bility  analysis  is  performed  (resonance  and  parametric 
instability).  To  solve  the  equation  for  random  load¬ 
ing,  the  Pierson-Moskowitz  power  spectrum,  describing 
the  wave  height,  is  first  transformed  into  an  approxi¬ 
mate  time  history  using  Borgman’s  method  with  slight 
modification.  The  equation  of  motion  is  then  solved, 
and  the  influence  on  the  tower  response  of  different  pa¬ 
rameter  values  such  as  buoyancy,  initial  conditions  and 
wave  height  and  frequency,  is  investigated. 


Key  words  :  Articulated,  Dynamics,  Random,  Sta¬ 
bility,  Chaos 


REVIEW  AND  PROBLEM  DEFINITION 

Compliant  platforms  such  as  articulated  towers  are  eco¬ 
nomically  attractive  for  deep  water  conditions  because 
of  their  reduced  structural  weight  compared  to  cmven- 
tional  platforms.  The  foundation  of  the  tower  does  not 
resist  lateral  forces  due  to  wind,  waves  and  currents;  in¬ 
stead,  restoring  moments  are  provided  by  a  large  buoy¬ 
ancy  force,  a  set  of  guylines  or  a  combination  of  both. 
These  structures  have  a  fundamental  frequency  well  be¬ 
low  the  wave  lower  frequency.  As  a  result  the  rela¬ 
tively  large  displacements,  geometric  nonlinearity  is  an 
important  consideration  in  the  analysis  of  such  a  struc¬ 
ture.  The  analysis  and  investigation  of  these  kind  of 
problems  can  be  divided  into  two  major  groups;  deter¬ 
ministic  and  random  wave  and/or  current  loading.  We 
briefly  review  wwk  in  this  area  in  the  next  two  subsec¬ 
tions. 

Deterministic  Loading 

Chakrabarti  and  Ck>tter  (1979)  [2]  analyzed  the  motion 
of  articulated  tower.  The  tower  is  articulated  by  a  uni¬ 
versal  joint  having  single  degree  of  freedom.  They  as¬ 
sumed,  linear  waves,  small  perturbations  about  an  equi¬ 
librium  position,  linear  drag  force  and  that  the  wind 
current  and  wave  are  colinear.  Their  resulting  equation 
of  motion  is, 

70  +  B(^)  +  Dii>  +  Ctl>  =  ( 1) 

where  I  is  the  total  moment  of  inertia  including  added 
mass,  B{i))  is  the  nonlinear  drag  term,  is  the  struc¬ 
tural  damping,  Cip  is  the  restoring  moment  due  to 
buoyancy  and  Mo  is  the  wave  moment.  An  analyti¬ 
cal  solution  is  then  compared  to  experimental  results 
to  show  good  agreement  as  long  as  the  system  is  inertia 
predominant,  and  not  drag  predominant. 


In  a  later  paper,  Chakrabarti  and  Cotter  (1980)  [3} 
investigated  transverse  motion,  the  motion  perpendic¬ 
ular  to  the  horu(Hital  velocity.  The  tower  pivot  is  as¬ 
sumed  to  have  two  angular  degrees  of  freedom  and  is 
taken  to  be  frictionleas.  It  is  also  assumed  that  the  mo¬ 
tion  is  not  coupled,  so  the  inline  solution  is  obtained 
(the  same  as  in  the  previous  paper)  from  which  the 
relative  velocity  between  the  tower  and  the  wave  is  ob¬ 
tained.  The  lift  f<»ce  (in  the  transverse  direction)  can 
then  be  obtained  and  the  linear  equation  of  motion  is 
solved  analytically  and  compared  to  experimental  re¬ 
sults.  The  comparison  shows  good  agreement,  espe¬ 
cially  when  the  drag  terms  are  small. 

Jain  and  Kirk  (1981)  [4]  investigated  the  dynamic 
response  of  a  double  articulated  offshore  structure  to 
waves  and  current  loading.  They  assumed  four  degrees 
of  freedom,  two  angular  for  each  Unk.  The  equations 
oi  motion  were  derived  using  Lagrange  equations.  In 
deriving  the  equations  of  motion  the  fdlowing  assump¬ 
tions  were  made:  drag  and  inertia  forces  tangent  to  the 
tower  are  negligible,  and  the  wave  and  current  veloci¬ 
ties  are  evaluated  at  the  upright  position  (small  angles 
assumption).  The  linearised  equations  were  solved  to 
find  the  natural  frequencies  of  the  system  and  then  nu¬ 
merically  solved  to  find  the  response  due  to  colinear  and 
mm-cdinear  current  and  wave  velocities.  They  found 
that  when  the  wave  and  the  current  velocities  are  co- 
linear,  the  response  of  the  top  is  sinusoidal,  while  for 
nonccdinear  velocities  the  response  is  a  complex  three 
dimensional  whirling  oscillation. 

Thompson  et  al.  (1984)  [5]  investigated  the  motions 
of  an  articulated  mooring  tower.  They  modeled  the 
structure  as  a  bilinear  oscillator  which  consists  of  two 
linear  oscillators  having  different  stiffnesses  for  each  half 
cycle, 

mx  +  CX  +  (ki ,  fcj)!  =  Fq  sin  ut,  (2) 

where  ki,  ate  the  stiffnesses  for  x  >  0  and  x  <  0 
respectively.  The  equaticm  is  solved  numerically  for 
different  spring  ratios  and,  as  expected,  harmonic  and 
subharnKHiic  resonances  appeared  in  the  respcmse.  A 
coiiq>ari8on  between  the  response  and  experimental  re¬ 
sults  of  a  reduced-scale  model  showed  good  agreement 
in  the  main  phenomenon. 

Choi  and  Lou  (1991)  [6]  have  investigated  the  be¬ 
haviour  of  an  articulated  offrhore  platform.  They  mod¬ 
eled  it  as  an  upright  pendulum  having  (me  degree  of 
freedom,  with  linear  springs  at  the  top  having  different 
stif&iesses  for  positive  and  negative  displacements  (bi¬ 
linear  oscillator).  The  equation  of  motion  is  simplified 
by  expanding  nonlinear  terms  into  a  power  series  and 
retaining  only  the  first  two  terms.  They  assumed  that 
the  ctxnbined  drag  and  inertia  ixunnent  is  a  harmonic 


function.  The  discontinuity  in  the  stiffness  is  assumed 
to  be  small,  and  thus  replaced  by  an  equivalent  contin¬ 
uous  function  using  a  least-squares  method  to  get  the 
following  DuflSng  equation 

10  +  c0 -i- ki$  +  kj$^ -i- ksff^  =  Mo  coBu)t,  (3) 

where  ki,  kj,  ks  are  spring  constants  depending  on 
buoyancy,  gravity  and  the  mooring  lines.  The  equation 
of  motion  is  solved  analytically  and  numerically,  and 
stability  analysts  is  performed.  The  analytical  solution 
agrees  very  well  with  the  numerical  solution.  The  main 
results  of  their  analysis  are  that  as  damping  decreases, 
jump  phenomena  and  higher  subharmonics  occur,  and 
chaotic  motion  occurs  only  for  large  waves  and  near  the 
first  subharmonic  (excitation  frequency  equals  twice  the 
fundamental  frequency);  the  system  is  very  sensitive  to 
initial  conditions. 

Seller  and  Niedzwecki  (1992)  [7]  investigated  the  re¬ 
sponse  of  a  multi-articulated  tower  in  planner  motion 
(upright  multi-pendulum)  to  account  for  the  tower  flex¬ 
ibility.  The  restoring  nooments  (buoyancy  and  gravity) 
were  taken  as  linear  rotational  springs  between  each 
link,  although  the  authors  say  that  nonlinear  springs 
are  more  adequate  for  this  model.  Each  Unk  is  assumed 
to  have  a  different  cross  section  and  density.  The  equa¬ 
tions  of  motion  are  derived  using  Lagrange’s  equations, 
in  which  the  generalized  coordinates  are  the  angular 
deflections  of  each  Unk.  The  equations  in  matrix  form 
we 

=  (4) 

where  [Af]  is  a  mass  matrix  that  includes  the  actual 
mass  of  the  Unk  and  added  mass  terms,  while  the  stiff¬ 
ness  matrix  [/if]  includes  buoyancy  and  gravity  effects. 
Damping  and  drag  forces  are  not  included  in  the  model. 
The  homogeneous  equations  for  a  tri-articulated  tower 
are  numerically  solved  to  study  the  effects  of  different 
parameters,  such  as  link  length,  material  density  and 
spring  stiffness,  on  the  natural  frequency  of  the  system. 

GotUb  et  al.  (1992)  [8]  analyzed  the  nonlinear  re¬ 
sponse  of  a  single  degree  of  freedom  articulated  tower. 
In  the  derivatim  of  the  equation,  the  expressions  for 
the  buoyancy  moment  arm,  added  mass  term,  and  drag 
and  inertia  moments  are  evaluated  along  the  station¬ 
ary  upright  tower  position  and  not  at  the  instantaneous 
position  of  the  tower.  The  governing  equation  is  of  the 
form 

0  +  y$  +  Rie)  =  Milt),  (5) 

where  /2(9)  =  a  sin  0  and  a  is  linear  function  of  buoy¬ 
ancy  and  gravity,  M{0,t)  is  the  drag  moment.  Approxi¬ 
mated  harmcmic  and  subharmonic  solutions  are  derived 
using  a  finite  Fourier  series  expansion,  and  stabiUty 


analysis  is  performed  by  a  Lyapunov  function  approach. 
The  solution  shows  a  jump  phenomenon  in  primary  and 
the  secondary  resonances. 

Random  Loading 

Muhuri  and  Gupta  (1983)  [9]  investigated  the  stochastic 
stability  of  a  buoyant  platform.  They  used  a  linear 
single  degree  of  freedom  model  as  follows 

i  +  2ci  +  (1  +  G(<))i  =  0,  (6) 

where  z  is  the  displacement,  c  is  the  damping  coeffi¬ 
cient  and  G(t)  is  a  stochastic  time-dependent  function 
due  to  buoyancy.  It  is  assumed  that  G(()  is  a  narrow- 
band  random  process  with  zero-mean.  A  criterion  for 
the  mean  square  stability  is  obtained  from  which  the 
!  lilowing  results  are  found:  for  c  >  1  the  system  is  al¬ 
ways  stable,  and  for  0  <  c  <  1  there  are  regions  of 
stability  and  instability. 

Datta  and  Jain  (1990)  [10]  investigated  the  response 

an  articulated  tower  to  random  wave  and  wind  forces. 
In  the  derivation  of  the  single  degree  of  freedom  equa¬ 
tion  of  motion  the  tower  is  dtscretized  into  n  elements 
having  appropriate  masses,  volumes  and  areas  lumped 
at  the  nodes,  and  there  is  viscous  damping.  The  equa¬ 
tion  of  motion  is, 

1(1  +  0(t))0  +  c0  +  Ril  +  vit))$  =  F(t),  (7) 

where  10{t)  is  the  time  varying  added  mass  term, 
Rv(t)  is  time  varying  buoyancy  moment  and  F(t)  is 
the  random  force  due  to  wave  and  wind.  The  Pierson- 
Moskowitz  spectrum  is  assumed  for  the  wave  height  and 
Davenport’s  spectrum  assumed  for  the  wind  velocity. 
The  equation  is  solved  in  the  frequency  domain  using 
an  iterative  method,  which  requires  that  the  deflection 
angle  $(i)  and  the  forcing  function  F(t)  be  decomposed 
into  Fourier  series:  The  coefficient  of  the  stn  and  cos 
are  then  found  iteratively.  From  their  parametric  study, 
they  concluded  the  following: 

1.  Ncmlinearities  such  as  large  displacements  and  drag 
force  do  not  influence  the  response  when  only  wind 
force  is  considered. 

2.  The  random  wind  forces  result  in  higher  responses 
than  do  wave  forces. 

3.  The  r  jn.8.  response  due  only  to  wind  forces  varies 
in  a  linear  fashion  with  the  mean  wind  velocity. 

In  a  later  paper,  Jain  and  Datta  (1991)  [11]  used  the 
same  equation  and  the  same  method  of  solution  to  in¬ 
vestigate  the  response  due  to  random  wave  and  current 


loading.  The  wave  loadings  (drag,  inertia  and  buoy¬ 
ancy)  are  evaluated  via  numerical  integration.  The  fol¬ 
lowing  results  were  obtained  from  the  parametric  study, 

1.  The  dynamic  response  is  very  small  since  its  fun¬ 
damental  frequency  is  well  below  the  wave’s  fun¬ 
damental  frequency. 

2.  Nonlinear  effects  (drag  force,  variable  buoyancy) 
have  considerable  influence  on  the  response. 

3.  Current  velocity  has  a  large  influence  on  the  re¬ 
sponse. 

Hanna  et  al.  (1983)  [12]  analyzed  the  nonlinear  dy¬ 
namics  of  a  guyed  tower  platform.  The  tower  is  rep¬ 
resented  by  a  lumped  parameter  model  consisting  of 
discrete  masses.  Each  mass  has  three  degrees  of  free¬ 
dom,  two  translations  and  one  rotation  about  the  ver¬ 
tical  axis.  The  external  forces  on  the  structure  are  ap¬ 
proximated  by  concentrated  forces  and  torques  at  the 
nodal  points.  The  equation  of  motion  is 

[Af]{ti}  +  [q{u)  +  [/!:(«)]{«}  =  {P(/.u,ii)},  (8) 

where  [Jl/]  is  the  total  mass  matrix  including  added 
mass  terms,  [C]  is  the  structural  damping  matrix  as¬ 
sumed  to  be  propcMTtional  to  the  mass  matrix  and  [/f(u)] 
is  the  total  nonlinear  stiffiiess  matrix  that  ir.cludes 
mooring  lines  effects,  soil  stiffness  and  geometric  stiff¬ 
ness.  {P(t,u,tt)}  is  the  nonlinear  dynamic  load  vec¬ 
tor  due  to  wave,  current  and  wind.  The  equation  is 
then  solved  using  normal  mode  superposition  and  the 
response  is  calculated  at  each  time  step.  This  method 
is  good  only  if  the  nonlinearities  are  not  large.  Deter¬ 
ministic  and  random  loading  are  considered.  The  so¬ 
lution  shows  insignificant  flexure  modes  while  the  tor- 
n<»al  one  has  a  noticeable  effect  on  the  deck  rotational 
response. 

Kanegaonkar  and  Haidar  (1988)  [13]  investigated  the 
nonlinear  random  vibration  of  a  guyed  tower.  They  in¬ 
cluded  nonlinearities  due  to  guylines  stiffnesses,  geome¬ 
try,  load  and  material.  The  simplified  planner  equation 
of  motion  is 

I0  +  e9  +  K9  +  kifi  =  M(t),  (9) 

where  A  is  a  spring  constant  depending  on  buoyancy, 
gravity  and  guyline  horizontal  stiffness,  and  ki  is  a  con¬ 
stant  depending  on  the  guyline  vertical  stiffness.  M{t) 
is  the  random  wave  loading.  The  equation  is  then 
solved  numerically  where  the  wave  height  is  defined  by 
the  Piersoa-Moskowitz  spectrum.  It  was  seen  that  the 
response  is  non-Gaussian  for  significant  wave  heights 
greater  than  5  m. 


Gerber  and  Engelbrecht  (1993)  investigated  the  re¬ 
sponse  of  an  articulated  mooring  tower  to  irregular  seas. 
It  is  an  extension  of  earlier  work  dcme  by  Thomps<ni  et 
al.  (1984)  [5].  The  tower  is  modeled  as  a  bilinear  oscil¬ 
lator,  that  is,  a  linear  oscUlator  with  different  stiffnesses 
for  positive  and  negative  deflections, 

mi-^ei  +  {ki,k2)x  =  F{t).  (10) 

Hie  random  forcing  function  F{t)  is  assumed  to  be  the 
sum  of  a  large  number  of  harmonic  c<«nponents,  each 
in  different  frequency,  a  procedure  similar  to  that  pro¬ 
posed  by  Borgman  (1969)  [14].  The  equation  is  then 
solved  analytically  since  it  is  linear  for  each  half  cycle. 
The  sdution  is  obtained  few  different  cases;  linear  oscil¬ 
lator  (both  stiffnesses  are  the  same),  bilinear  oscillator, 
and  for  the  case  of  impact  oscillator  (a  rigid  cable)  in 
which  oscillation  can  occur  mly  in  one  half  of  the  cy¬ 
cle.  For  future  study  they  suggest  to  include  nonlinear 
stiJfiiess  and/or  using  a  different  spectrum  to  describe 
the  wave  height. 

Problem  Definition 


PROBLEM  DESCRIPTION 

A  schematic  of  the  structure  is  shown  in  Fig.  1.  It 
consists  of  a  tower  submerged  in  a  fluid  having  a  con¬ 
centrated  mass  at  the  top  and  one  degree  of  freedom 
0  about  the  z  axis.  The  tower  is  subjected  to  wave 
loading.  Two  coordinate  systems  are  used;  one  fixed 
(x,y,z)  and  the  second  attached  to  the  tower  (x'.y'.r'). 
AU  forces/moment  are  derived  in  the  fixed  coordinate 
system,  which  means  that  the  tower  rectilinear  velocity 
is  resolved  into  x,  y  coordinates.  The  motion  td  the 
tower  is  assumed  to  be  only  in  plane. 

This  problem  has  similarities  to  that  of  an  inverted 
pendulum,  but  due  to  the  presence  of  gravity  waves, 
additional  considerations  are  made  : 

1.  Forces  due  to  buoyancy  and  vertical  wave  velocity 
are  summed  wd  denoted  as  Tq. 

2.  Drag  fcxrces  proportional  to  the  square  of  the  relar 
tive  velocity  between  the  fluid  and  the  tower  need 
to  be  considered. 


In  this  paper,  the  respemse  of  an  articulated  tower  sub¬ 
merged  in  the  ocean  is  investigated.  The  nonlinear  dif¬ 
ferential  equation  of  motion  is  derived,  including  nonlin¬ 
earities  due  to  gecxnetry,  coulomb  damping,  drag  force, 
added  mass  buoyancy.  All  forces/moments  are  eval¬ 
uated  analytically  and  explicitly  at  the  instantaneous 
position  of  the  tower  and,  therefore,  they  are  time- 
dependent  and  highly  nonlinear.  The  equatiem  is  then 
numerically  solved  using  ‘ACSL’  •  Advanced  Continu¬ 
ous  Simulation  Language  [15],  a  powerful  sefftware  lan¬ 
guage  for  deterministic  and  random  wave  loading  us¬ 
ing  the  Pierson-Moskowitz  wave  height  spectrum.  A 
harmonic  and  subharmonic  solutions  for  deterministic 
wave  heights  are  obtained.  The  response  to  random 
wave  heights  for  different  significant  wave  heights  is  in¬ 
vestigated,  the  influence  of  coulomb  damping  on  the 
response  is  analyzed,  and  chaotic  regimes  of  behavior 
are  identified. 

The  distinctiems  between  this  study  and  the  literature 
with  which  we  are  aware  are  that, 

1.  A  sound  and  exact  derivation  of  the  nonlinear 
equation  of  motion  is  provided. 

2.  All  terms  in  the  equation  of  motion  are  analytically 
derived. 

3.  Coulomb  friction  in  the  tower  hinge  is  added. 

4.  Usage  ‘ACSL’  for  the  numerical  solution  pro¬ 
vides  an  easy  way  to  modify  parameters  and  per¬ 
form  a  sensitivity  study. 


3.  Fluid  inertia  forces  due  to  fluid  acceleration  are 
part  of  the  loading  environment. 

4.  Fluid  added  mass  is  directly  included  in  the  inertia 
forces. 


EQUATION  OF  MOTION 


The  equation  of  motion  is  derived  using  Newton’s  sec¬ 
ond  law,  setting  all  moments  (external  and  inertial)  act¬ 
ing  on  the  tower  about  hinge  to  sero.  We  will  find  the 
equation  of  motion  to  be  of  the  form, 

mm = (t,w./(fl),fl,siyn(fl)) ,  (11) 

where,  is  the  effective  position  dependent  mo¬ 

ment  of  inertia,  ff(ff),  f(9)  are  nonlinear  functions  of  9 
(trigcmometric  functions),  u  is  the  wave  frequency  and 
M  is  the  sum  of  all  external  moments  that  act  on  the 
tower.  Certain  assumptions  have  been  made  in  deriv¬ 
ing  the  nonlinear  equation  of  motion.  These  are  listed 
below. 

Assumptions 

1.  The  tower  stiflhess  is  infinite:  El  =  oo. 

2.  The  hinge  has  coulomb  damping. 

3.  The  tower  has  a  uniform  mass  per  unit  length,  m 
and  is  of  length  /  and  diameter  D. 

4.  The  tower  is  a  smooth  slender  structure  with  uni¬ 
form  cross  section. 

5.  The  end  mass  M  is  considered  to  be  concentrated 
at  the  end  of  the  tower.  (It  has  no  volume.) 

6.  The  tower  length  is  greater  than  the  fluid  depth, 
but  the  dynamics  is  not  limited  to  the  case  of  Af 
always  being  above  the  mean  water  level. 

7.  The  structure  is  statically  stable  due  to  the  buoy¬ 
ancy  force. 

8.  The  waves  are  linear  having  random  height. 

9.  Morison’s  fluid  force  coefiBcients  Co  and  Cm  are 
constant. 

10.  The  center  of  mass  (c.g.)  of  the  tower  is  at  its 
gecmietric  center. 

11.  Currents,  wind  and  wave  slamming  forces  are  not 
included. 

Farces/Moments  Acting  on  the  Tower 

Fig.  1  describes  the  external  forces  acting  on  the  tower. 
These  are  : 


2.  Fv,  Fh  are  the  vertical  and  horizontal  fluid  forces 
due  to  fluid  drag  and  inertia. 

3.  Mg,  mlg  are  the  forces  due  gravity. 

We  next  describe  these  forces  and  moments. 

Inertia  Moment 

The  inertia  moment  equals  the  total  moment  of  inertia 
of  the  tower  plus  the  fluid  added  mass  term,  multipUed 
by  the  angular  acceleration  of  the  tower, 

Mt  =  (Jo  +  MP)9  +  Mj,,  (12) 

where  Jo  is  the  moments  of  inertia  of  the  tower  about 
point  ‘o’  and  M//  is  the  fluid  added  mass  moment.  As¬ 
suming  that  the  fluid  added  mass  due  to  the  end  mass 
is  negligible,  the  total  inertia  moment  is 

Mj  =  Jeffl  (13) 

where  Jtj j ,  the  total  moment  of  inertia,  is 

Jli  =  liml  +  M)l^+  (14) 

~CaP^^L\1  +  tan^9). 

where  L  is  the  projecti(»i,  in  the  x  direction,  of  the 
submerged  part  of  the  tower  and  is  quantified  later.  Ca 
is  the  added  mass  coefficient  which  equals  Ca  =  Cm~^> 
where  Cm  is  the  inertia  coefficient. 

Moments  due  to  Gravity,  Buoyancy 

The  fdlowing  moment  b  due  to  terms  that  do  not  de¬ 
pend  on  the  fluid  velocity,  such  as  gravity  and  buoyancy 
forces, 

Mf  =  Tolt  —  {M +  m^)glsin9,  (15) 

where  g  is  the  gravitational  constant.  To  is  the  buoyancy 
force  which  is  time  dependent  and  /»  is  its  moment  arm; 

D* 

To  =  pgx—L,,  (16) 

4 

where  Lg,  which  is  the  length  of  the  submerged  paer  of 
the  tower,  is 

where  d  is  the  mean  water  level,  p  is  the  fluid  den¬ 
sity  is  the  wave  height  elevation  evaluated  at  the 
mean  water  level, 

ri{9,t)  =  ~Hcos{kdtaxi9  —  ut  +  e),  (18) 


1.  To  is  a  vertical  buoyancy  force. 


H  is  the  wave  hei^t,  u  is  the  wave  frequency,  k  is  the 
wave  number.  The  buoyant  force  acts  at  the  center  of 
mass  of  the  submerged  part  of  the  tower  U, 

=  j^tan’tfcostf+ii, +  ^^tan’tf8mtf  (19) 

Morison’s  Equation  for  Wave  Force 

In  general,  the  fluid  forces  acting  on  a  slender  tower 
are  of  two  types:  drag  and  inertia.  The  drag  force 
is  proportional  to  the  square  of  the  relative  velocity 
between  the  fluid  and  the  tower,  and  the  inertia  force 
is  proportimial  to  the  fluid  acceleration.  The  drag  and 
inertia  forces  per  unit  length  are  given  by  Morison’s 
equation, 


When  solving  the  equation  of  motion  numerically,  the 
sign  of  the  relative  velocities  and  of  (d  -  IcobO)  are 
checked  in  each  time  interval.  This  is  done  to  set  the 
correct  direction  of  the  drag  moments  and  the  limits  of 
integration,  respectively. 

Damping  Moment 

The  tower  hinge  is  assumed  to  dissipate  energy  via 
coulumb  friction.  In  this  section  the  damping  moment 
is  evaluated.  The  general  form  of  the  friction  fwce  is 

Ffr  =  N  ii[8gn($)],  (25) 

where,  fi  is  the  friction  coefficient,  and  N,  the  normal 
force,  is 


r/i  =  c*,|u-v|(u-v)  +  c„„u,  (20) 

where, 


Cdd  = 

Cmm  —  CiipV—',  (21) 

Fyt  is  the  fluid  force  per  unit  length,  V  is  the  tower 
absdute  velocity  and  U  is  the  fluid  absolute  velocity 
evaluated  at  the  instantanuouse  position  of  the  tower. 
The  moments  due  to  the  fluid  velocity  and  acceleration 
are  as  follows.  The  inertia  moments  are 


Mn  = 

Mj,  = 


itxdx 


wx  tan  fidx. 


(22) 


and  the  drag  moments, 


Afj>h  =  Cid  f  !«-«,!(«-  v^)xdx  (23) 
Jo 

Mdv  =  Cdd  I  I  u;  — I  (tu  -  t>,)ztantfdi. 

Jo 


where  A//a,  Jtf/v  are  the  inertia  moments  due  to  hori¬ 
zontal  and  vertical  accelerations,  respectively  and  Mdh, 
Mdv  ak  the  drag  moments  due  to  horizontal  and  ver¬ 
tical  relative  velocities,  respectively.  L,  which  is  the 
upper  limit  of  the  integral,  depends  cm  the  angle  0  as 
follows  : 


f  lcoB$  ifd>/cos0 
\  d  if  d  <  lco6$. 


(24) 


This  means  that  the  moments  are  calculated  only  for 
the  submerged  part  of  the  tower.  The  integrals  in  equar 
tknis  (22)  and  (23)  are  evaluated  analytically  using  the 
aymbolic  software  ‘MAPLE’. 


JV  =  cos  -1-53  F,  sin  fl,  (26) 

where  ^  Fn  forces  in  the  z,  y  di- 

recti(His,  respectively.  These  forces  are  due  to  gravity, 
buoyancy  and  fluid  drag  and  inertia.  If  we  assume  a 
hinge  radius  Rk,  then  the  damping  moment  is 

V/.  =  {Ef-  coefl-l-^F,  sintf^  • 

RkP[B9nm  (27) 

Finally,  adding  all  the  external  moments,  equations 
(15),  (22),  (23),  (27),  and  setting  them  equal  to  the  to¬ 
tal  inertia  moment,  equation  (13),  yields  the  governing 
nonlinear  differential  equaticm  of  motion, 

JeJJ0=  ~Mjr  —  Mf  +  Mlh  —  Mjv  +  Moh  —  Mdv  (28) 

Both  sides  of  this  equation  are  a  nonlinear  functions  of 
9,  9,  t  and  u. 


NUMERICAL  SOLUTION 


The  numerical  solution  for  the  nonlinear  differential 
equation  (28)  was  performed  using  ‘ACSL’  and  the 
analyses  of  the  results  was  performed  using  ‘MATLAB’. 
The  ‘ACSL’  code  written  for  this  application  is  avail¬ 
able  upon  request.  In  order  to  build  confidence  in  the 
solution,  some  test  cases  have  been  solved  using  the 
following  physical  parameters  : 


Tower  properties 


1.  /  -  Length  of  the  tower  =  400  (m) 


2.  D-  Tower  diameter  =  15  (m) 

Z.  M  -  End  mass  =  2.5  *  10*^  (Kg) 

4.  m  •  Tower  mass  per  unit  length  =  2  *  10^  (^) 


5.  /I  •  FMctioo  coefficient  =  0.1-0.4 

6.  Aft  •  Hinge  radius  =  1.5  (m) 

Fluid  properties 

1.  d  -  Mean  water  level  =  350  (m) 

2.  Cd  -  Drag  coefficient  =0.6  •  1.0 

3.  Cm  '  Inertia  Coefficient  =1.5 

4.  p  -  Water  density  =  1025  (^) 

5.  u  -  Wave  frequency  =  0.2  -  1.0  (^) 

Respoiise  for  Detenninistic  Wave  Heights 

In  this  section  the  response  of  the  tower  to  detenninistic 
wave  heights  is  established.  Free  vibration  and  stability 
analyses  are  performed. 

First,  the  natural  frequency  of  the  tower  is  found  as 
if  it  was  an  upright  pendulum  subjected  to  a  constant 
tension  force  To  and  gravitational  force  Mg.  Fig.  2 
shows  the  tower  response  to  an  impulse  for  Cn  =  0  and 
A  =  0  in  the  tiiiK  and  frequency  domain.  The  response 
is  harmonic  in  the  tower  natural  frequency  n„  =  0.026 
(Hx),  and  it  agrees  with  the  calculation  for  the  natural 
frequency  of  a  pendulum  : 

=  (29) 

Fig  3  (a,b)  shows  the  same  but  for  p  =  0  and  Cd  = 
1,  the  decay  here  is  not  linear  since  the  drag  fmce  is 
propotional  to  the  velocity  squared.  Fig  3  (c,d)  shows 
the  free  vibration  with  frictional  damping,  /i  =  0.1,  in 
the  time  and  frequency  domain.  Here  the  amplitude 
decays  linearly  with  time  as  expected  when  coulomb 
damping  is  present.  From  the  figure  we  can  evaluate 
the  equivalent  damping  ratio  f<»  Cj>  =  0  to  be  (  =  0.02. 
Since  the  damping  is  nonlinear,  the  ‘natural  frequency’ 
and  its  multipliers  are  seen  in  the  figure. 

Fig.  4  shows  the  response  in  the  time  domain  and 
frequency  dcnnain  for  A  =  1  (m)  and  wave  loading  fre¬ 
quency  u  =  0.064  {Hz).  Figures  (a,  b)  are  with  p  =  0 
and  (c,  d)  are  with  ft  =  0.4.  In  the  frequency  domain, 
the  tower  ‘natural  frequency’  and  the  wave  natural  fre¬ 
quency  are  clearly  seen.  It  can  be  seen  that  the  friction 
has  a  damping  and  stabilizing  effect  on  the  system. 
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Figure  2;  Tower  Natural  Frequency,  Cd  ~  0 
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Figure  3:  Tower  FVee  Vibration  with  Damping,  (a,  b) 
ft  =  0.1,  (c,  d)  Cd  =  1 
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FigUK  4:  Tower  response  to  wave  excitation,  Time  and 
Frequency  Domain,  H  n  I  (m) 


Since  this  is  a  nonlinear  system,  its  response  to  har¬ 
monic  excitation  at  the  system’s  ‘natural  frequency’  or 
its  multipliers  can  be  multivalued,  indicating  the  oc¬ 
currence  of  a  ‘jump’.  Fig.  5  shows  the  tower  respcxise 
to  harimmic  wave  excitation  at  w  =  (2„  (harmonic  so¬ 
lution)  with  and  without  damping.  If  the  system  was 
linear,  such  loading  would  have  caused  resonance  instar 
bility.  But  since  the  system  is  nonlinear,  an  amplitude 
change  (beating)  can  be  seen  Fig.  5  (a),  indicating  that 
‘jumps’  occur  from  one  amplitude  to  the  other  (Wilson 
[16]  pp.  147).  It  can  also  be  seen  *  .it  although  the  ex¬ 
citation  frequency  is  constant,  the  frequency  response 
is  not  Fig.  5  (b).  Again  the  reason  is  the  nonlinear 
system  characteristics.  FVom  the  response  with  friction 
(figures  (c,  d))  we  see  that  the  response  is  not  smaller 
since  it  is  unstable,  but  the  ratio  between  the  large  and 
small  amplitude  of  the  beating  gets  smaller. 

Fig.  6  shows  the  response  due  to  harmonic  wave  ex¬ 
citation  at  twice  the  system  ‘natural  frequency’  (sub¬ 
harmonic  solution),  with  (c,  d)  and  without  friction  (a, 
b).  We  can  see  the  beating  phenomenon,  and  from  the 
frequency  domain  we  see  that  the  response  consists  of 
the  ‘natural  frequency’and  it  multipliers.  This  is  due  to 
the  nonlinear  nature  of  the  system. 


Figure  5;  Tower  response  to  harmonic  wave  excitation 
at  the  ‘Natural  Frequency’  -  Beating  Phenomenon 
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Figure  6:  Tower  response  to  sub-harmonic  wave  excita¬ 
tion  at  twice  ‘Natural  Frequency’ 


The  direction  of  the  drag  force  is  opposite  to  the  di¬ 
rection  of  the  relative  velocity  between  the  tower  and 
the  wave.  If  the  direction  of  the  relative  velocity  is  the 
same  as  the  tower  velocity,  the  drag  will  always  stabi¬ 
lise  the  system.  The  reason  is  that  the  drag  moment 
will  always  be  opposite  to  the  direction  of  the  deflec¬ 
tion  angle  B.  But  if  the  relative  velocity  direction  is 
(^posite  to  tower  velocity,  the  drag  moment  can  cause 
iiMtability.  For  aero  wave  velocity  the  drag  force  always 
stabilises  the  system,  as  can  be  seen  from  Fig.  3.  Fig. 
7  describes  the  horisontal  and  vertical  components  of 
the  relative  velocity  and  tower  velocity. 
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Figure  7:  Relative  and  Tower  Velocities  -  (a,  c)  Hori¬ 
sontal  and  (b,  d)  Vertical 


Chaotic  Response 

Because  of  the  nrnilinear  characteristics  of  the  system, 
the  governing  differential  equation  of  motion,  chaotic 
motkm  can  occure  under  certain  conditions.  Fig.  8 
shows  the  response  of  the  tower  to  a  harmonic  wave 
excitation  at  a  frequency  u  =  O.flSRn  •  We  see  that  the 
response  'jumps’  between  two  stable  equilibrium  states. 
This  phenomenon  can  imply  a  chaotic  system. 
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Figure  8:  Tower  response  to  harmonic  wave  excitation 
at  the  ‘Natural  Frequency’,  Co  =  0.6 


Moon  [17]  provides  several  ways  to  identify  chaotic 
vibrations.  We  next  use  three  ways  to  prove  that  the 
system  under  consideration  is  chaotic.  Fig.  9  describes 
the  response  of  the  tower  to  a  harmonic  wave  excitation 
having  a  frequency  ofcv  =:  0.09  {Hz).  Fig  9  (a,  b)  shows 
the  response  in  the  time  and  frequency  domains  and  the 
phase-plane  trajectory  (c).  From  the  time  history  it  is 
clearly  seen  that  the  system  jumps  between  two  stable 
equilibrium  states.  The  frequency  domain  shows  mul- 
tihammnic  energy  although  the  system  has  one  degree 
of  freedom;  another  indicator  of  chaos.  Finally,  from 
the  phase-plane  trajectory  we  see  that  the  orbits  never 
close  or  repeat.  Thus,  we  conclude  that  this  model  has 
the  characteristics  of  a  chaotic  system.  This  chaotic 
phenomenon  appeared  only  when  the  wave  height  was 
larger  than  10  (m). 
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Figure  13:  Tower  Response  •  (a,  b)  Deflection  angle 
Figure  11:  Tower  Response  -  Deflection  angle  and  Total  9{t  =  0)  =  0.05  (rad)  and  (c,  d)  ff(i  =  0}  =  0.01 
Drag  Moment  tor  If t  —  9  (m)  and  Co  =  0.6  (dashed  (rad/aec). 
line)  and  1  (  solid  line). 


DISCUSSION  AND  SUMMARY 
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Figure  12:  Tower  Response  -  Deflection  angle  -  Time 
and  Ftequency  Domain  for  Ht  —  9  (m)  and  Co  =  1> 
(a,  b)  p  =  0  and  (c,  d)  it  =  0.4 


The  nonlinear  differential  equation  of  motion  for  an  ar¬ 
ticulated  tower  submerged  in  the  ocean  is  derived.  Geo¬ 
metric  as  well  as  force  nonlinearities  are  included  in  the 
derivation.  The  wave  velocities  and  accelerations  are 
determined  at  the  instantaneous  position  ci  the  tower, 
a  fact  that  added  to  the  nonlinearities  of  the  equation. 
The  equation  is  sdved  numerically  using  ‘ACSL’  for 
deterministic  and  random  wave  loading. 

The  response  of  the  tower  to  harmonic  wave  excita¬ 
tion  at  its  ‘natural  frequency’  and  at  twice  its  ‘natural 
frequency’  demonstrates  beating,  where  the  amplitude 
varies  between  two  extremes.  This  beating  is  due  to 
the  ncmlinear  behaviour  of  the  system.  Coulomb  damp¬ 
ing  reduces  the  beating  phenomenon  and  the  response 
amplitude,  so  it  has  a  stabilizing  effect  on  the  system. 
When  the  system  is  excited  at  an  arbitrary  frequency, 
and  the  wave  height  is  greater  than  about  10  (m),  the 
response  ‘jumps’  between  two  stable  equilibria,  exhibit¬ 
ing  chaotic  behaviour. 

To  solve  the  equation  for  random  wave  loading,  the 
Pierson-Moskowitz  spectrum  that  describes  the  wave 
height  distribution  was  first  transformed  into  a  time 
history.  The  equation  was  solved  for  two  significant 
wave  heights.  Again  the  response  was  periodic  consist¬ 
ing  of  the  fundamental  tower  frequency  and  its  multipli- 
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Figure  9:  Tower  respcmse  to  a  harmonic  wave  excita¬ 
tion,  u  =  0.09  (Hz)  -  Chaos 


0.2  (rad/see)  articulated  tower  are  designed  to  have  a 
‘natural  frequency’  lower  than  that. 
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Random  Wave  Height 

In  this  section  the  tower  response  to  random  wave 
height  excitation  is  investigated.  The  wave  height  dis¬ 
tribution  is  generally  expressed  in  the  form  of  a  power 
spectral  density.  For  a  simulation  of  the  respcxise 
in  the  time  dmnain,  the  Pierson-Moskowitz  spectrum, 
which  is  assumed  to  describe  the  spectrum  cf  the  wave 
height,  is  transformed  into  a  an  infinite  series  of  time- 
dependent  harmonic  functi<»i8.  This  is  accomplished 
using  a  method  by  B<»gman  [14],  described  also  in  Wil¬ 
son’s  book  [16]. 

Response  for  Random  Wave  Height 

In  this  section,  the  influence  of  different  significant  wave 
heists,  drag  coefficient  ,  hinge  fricti<m  and  non-zero 
intial  condition  on  the  response  is  investigated. 

Fig.  10  compares  the  tower  response  for  H,  ^  9 
and  15  (m).  Here  we  enlarged  the  buoyancy  force  so 
that  the  tower  ‘natural  frequency’  is  On  =  0.07  {Hz) 
=  0.44  {rad/aec).  It  can  be  seen  that  the  deflection 
angle  fw  if,  =  9  (m)  is  of  the  same  magnitude  as  for 
for  Ht  =  15  (m),  although  the  later  maximum  wave 
height  is  four  times  larger  then  the  fonner.  The  reason 
is  that  the  natural  frequency  of  the  tower  coinsides  with 
where  most  of  the  spectral  energy  for  if,  =  9  (m)  is 
located.  Since  the  lowest  wave  freqency  is  about  u  sz 


Figure  10:  Tower  Response  •  Deflection  angle  and  Total 
Moment  for  H,  ss  15  (dashed  line)  and  9  (solid  line)  (m) 

A  comparison  of  the  tower  response  for  different  val¬ 
ues  of  the  drag  coefficient  are  plotted  in  Fig.  11.  This 
figure  shows  the  deflection  angle  and  the  total  drag  mo¬ 
ment,  Md  —  Mph  —  Afpvi  for  Cp  =  0.6  and  1,  and 
H,  =  9  (m).  It  can  be  seen  from  the  figure  that  the 
deflection  angle  and  the  drag  moment  are  larger  for 
Cp  =  1,  since  the  drag  moment  in  the  Morison  equa¬ 
tion  is  proportional  to  the  drag  coefficient.  The  results 
were  similar  fw  different  significant  wave  heights. 

The  friction  effect  on  the  response  is  shown  clearly  in 
Fig.  12.  Here  the  significant  wave  height  is  H,  =  9  and 
Cp  =  1.  It  is  seen  that  the  response  is  much  smaller  and 
smocAher  due  to  the  additional  friction  in  the  system. 

The  response  for  nonzero  initial  conditions  is  depicted 
in  Fig.  13.  The  response  for  9(t  =  0)  =  0.01(rad/sec) 
(a,  b),  and  9(i  =  0)  =  0.05(rad)  (c,  d).  The  drag 
coefficient  is  Cp  =  0.6  and  H,  =9  (m). 


era.  For  significant  wave  height  of  9  (m),  the  response 
was  larger  than  that  for  15  (m),  since  in  the  former  the 
tower  ‘natural  frequency’  coincides  with  the  frequen¬ 
cies  where  most  of  the  energy  is  located.  The  response 
with  coulomb  damping  shows  that  friction  stabilises  the 
system.  Notice  that  in  order  to  reduce  stresses  in  the 
structure,  the  friction  moment  has  to  be  low  enough  so 
that  the  tower  can  comply  with  the  wave  loading. 

A  mwe  realistic  model  having  two  angular  degrees  of 
freedom  is  being  analyzed  at  the  present  time.  The  re¬ 
sponse  due  to  wave,  current  (colinear  and  ncm-colinear) 
and  vortex  shedding  loading  is  investigated  loading  and 
results  will  be  published  in  the  near  future. 
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Abstract 


In  Rodriguez  and  Van  Kampen’s  1976  paper  [1],  a  method  of  extracting  information  from 
the  Fokker-Planck  equation  without  having  to  solve  the  equation  is  outlined.  The  Fokker- 
Planck  equation  for  a  Duffing  oscillator  excited  by  white  noise  is  expanded  about  the  intensity 
of  the  forcing  function,  a.  This  expansion  is  carried  to  order  However  no  studies 

are  made  of  the  effects  of  the  order  of  the  expansion,  variation  of  the  parameters,  nor  are 
comparisons  made  to  experimental  results.  In  this  paper,  the  expansion  is  carried  to  a  higher 
order,  C7(a§),  results  are  presented  and  compared  to  Monte- Carlo  experiments  using  both 
white  and  colored  noise,  and  parametric  studies  are  performed  on  the  intensity  of  the  forcing 
function  and  the  damping  coefficient.  It  is  found  that  the  expamsion  method  works  well  for  the 
case  of  white  noise  and  for  colored  noise  where  the  correlation  time  is  less  than  0.1  seconds, 
but  falls  to  give  certain  detadls.  It  is  also  found  that  the  system  behaves  as  expected  when  the 
parameters  are  varied. 
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1  Introduction 


The  Fokker>Planck  equation  has  proven  to  be  a  useful  tool  in  the  analysis  of  simple  nonlinear 
oscillators  excited  by  stochastic  processes.  As  a  partial  differential  equation  for  the  probability 
density  function  of  the  response,  its  solution  completely  defines  the  solution  of  the  problem. 
It  can  be  used  to  analyze  both  a  single  oscillator  of  the  form 

mi  +  7(i,  i)i  +  A:(z,  x)x  =  .F(t),  (1) 

or  a  system  of  multiple,  linked,  oscillators  of  the  form 

Mi  +  r(x,x)x  +  K(x,x)i  =  Z(t).  (2) 

In  many  cases,  a  physical  system  can  be  approximated  by  such  a  system  of  nonlinear 
oscillators.  The  systems  so  modeled  can  range  from  a  Brownian  particle  to  structures  excited 
by  von  Karman  vortex  shedding.  Such  modeling  can  be  useful  for  gaining  insight  into  a 
problem  and  the  way  in  which  the  system  will  behave  as  certadn  parameters  are  varied. 

Once  one  has  decided  on  the  system  of  oscillators  to  be  used  to  represent  the  physical 
system,  the  derivation  of  the  Fokker-Planck  equation  is  relatively  straightforward,  although 
tedious.  The  problem  remains  of  how  to  solve  it  for  the  probability  distribution  of  the  response. 
In  a  very  few  cases,  the  Fokker-Planck  equation  can  be  solved  analytically,  but  in  most  cases 
no  analytical  solution  exists  and  one  usually  must  resort  to  a  numerical  solution.  However  this 
can  be  computationally  intensive  and  ^ves  little  insight  into  the  larger  problem. 

In  their  1976  paper,  Rodriguez  and  van  Kampen  outline  a  method  of  dealing  with  the  case 
of  an  oscillator  excited  by  weak  Gaussian  white  noise.  The  Fokker-Planck  equation  of  the 
system  is  expanded  about  the  intensity,  a,  of  the  driving  function.  This  expansion  is  carried 
to  the  order  0(aa ).  In  this  way  the  statistics  of  the  fluctuations  cam  be  obtained  directly.  This 
method  shows  promise  as  a  way  to  use  the  Fokker-Planck  equation  to  gain  useful  information 
about  a  wider  variety  of  systems  than  was  p<wsible  before. 

This  is  the  first  of  a  planned  series  of  papers  exploring  the  usefullness  of  this  method.  As 
in  the  original  paper,  the  method  is  applied  to  the  problem  of  a  Duffing  oscillator  excited  by 
Gaussian  white  noise.  The  inherent  assumptions  of  the  method  are  explained  here  in  detail. 
The  expansion  is  carried  both  to  the  same  order  as  in  the  original  paper,  and  to  order  0{a3 ). 
Results  are  presented  and  compared  to  a  Monte  Carlo  experiment.  Parametric  studies  were 
done  on  the  parameter  of  expansion  as  well  as  on  the  other  important  variable  in  the  expansion: 
the  coefficient  of  damping. 

2  Expansion  of  the  Fokker  Planck  Equation  for  a 
Duffing  Oscillator 

As  in  the  Rodriguez  and  van  Kampen  paper,  heretofore  called  the  “Original  paper,”  the 
system  under  consideration  is  a  Duffing  oscillator  in  a  heat  bath.  The  equation  of  motion  can 
be  written  simply  as 

X  -1-  7X  -I-  X  -I-  x^  =  F{t).  (3) 

F{t)  is  a  Langevin  force  [2]  and  is  assumed  to  be  Gaussian  white  noise  with  the  following 
properties: 

<Fit)>  =  0 

<F{t)F{t')>  =  2a6it-t').  (4) 
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It  is  assumed  that  the  immersion  of  the  oscillator  takes  place  at  time  t  —  0  and  that  the  system 
is  not  necessarily  at  rest. 

f{x,  V',  t)  is  defined  as  the  joint  probability  density  function  of  x  and  v  =  x  at  time  t.  This 
leauls  to  the  following  Fokker-Planck  equation  (see  Ochi  [3]  for  a  complete  derivation): 


«/  .  Jf 


d 


d^f 


at  '  dx 


dv^' 


(5) 


As  /(z,  v;  t)  is  a  complete  description  of  the  system  response,  solution  of  the  above  differential 
equation  for  /  constitutes  a  solution  of  the  problem. 

The  same  assumptions  on  the  sizes  of  the  variables  in  Eq.  5  are  made  as  in  the  original 
paper:  namely  that  7,z  and  v  will  all  be  assumed  to  be  of  the  same  order  of  magnitude,  and 
much  larger  than  a. 

In  the  original  paper,  it  is  assumed  that  the  system  response  due  to  the  forcing  function  will 
be  small  as  compared  to  the  deterministic  response  due  to  the  initial  conditions.  Therefore, 
the  total  response  can  be  viewed  as  random  fluctuations,  and  Av,  superimposed  onto  the 
deterministic  response  to  the  initial  conditions.  Furthermore,  the  random  fluctuations  will 
be  of  order  a^.  Because  the  only  source  of  energy  is  F{t),  the  power  input  to  the  system  is 
proportional  to  a.  But  if  the  system  is  to  remain  stable,  then  the  viscous  power  dissipation 
of  the  fluctuations  caused  by  the  influence  of  F{i)  must  be  of  equal  average  magnitude  as 
the  power  input.  Therefore  Oi'f  rj^)  =  0(a)  or  0{ti)  =  0(aa).  But  the  kinetic  energy  of 
the  fluctuations  must  be  of  the  same  order  as  the  potential  energy.  So  O(C^)  =  0{ri^)  which 
implies  that  0(C)  =s  0(ai).  Therefore,  A*  =  or^C  aJid  Av  =  Of^»7  where  C  and  t)  are  of  order 
unity.  Therefore  the  following  substitutions  are  made: 

X  =  ^t)  +  aaC.  (6) 

V  =  +  (7) 


In  the  original  paper,  the  initial  conditions,  ^0),  and  ^^(0)*  are  assumed  to  be  zr  )  and 
the  expansion  carried  through  to  give  the  time  derivatives  of  the  second  order  moments.  In 
Weinstein  and  Benaroya,  [4],  the  expansion  is  explained  in  somewhat  greater  detail  and  the 
fourth  order  moments  are  a^  derived.  In  total  the  following  eight  equations  are  derived: 

^  <  C^  >=  2'<  C»7  >  (8) 

>=<  r]^>  -  >  -1  <CV>  >  +  C?(at)  (9) 

4:  <  >=  -2  <  (r]  >  -27  <  >  -2q  <  C^»7  >  +2  +  C7(q*  )  (10) 

±<c'>=4<e’,>  (11) 

|<(’|)>=3<CV>  -  <C\  >-<<">  + 0(a>)  (12) 

4:  <  (W  >=  2  <  Cl®  >  -27  <  iW  >  -2  <  C’l?  >  +2  <  >  +  O(oi)  (13) 

4  <  C»7®  >=<  V*  >  -37  <  >  -3  <  CW  >  +6  <  c»?  >  +  0{a^ )  (14) 

at 

4  <rt*>=  -47  <T)*>  -A  <  C»?®  >  +12  <  7®  >  +  0(at ).  (15) 

at 
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If  it  is  the  stationary  behavior  that  is  of  interest,  then  one  can  obtain  the  equilibrium 
values  of  these  quantities  by  setting  each  time  derivative  equal  to  zero.  Then  the  equilibrium, 
or  stationary,  states  of  each  expectation  can  be  found  by  simple  linear  algebra: 


<C^>e,  =  ^-0^ 

<C®»?>e,  =  0+  C7(Qt) 

<C»7>t,  =  0 

<(,W>tq  =  :^+  0{a\) 

<»7^>e,  =  A+  C?(q3) 

<CT/3>e,  =  0+  C?(ot) 

<C^>.,  =  :^+  O(o^) 

<n*>tq  =  :^+  C?(q*). 

If  one  is  interested  in  the  transient  response  of  the  system,  and  if  one  can  accept  a  solution 
of  order  C7(a^),  then  one  can  obtain  it  an^yticaily.  Elquations  8  to  10  can  be  written  in  matrix 
form  as  follows: 


dt 


«  C" »  ^ 

<<  C*? »  =  M3 

«  ri^  »  / 


I  0  2 

where  M3  =  <  -1  —7 
0  -2 


«  c* » 
«  C»7» 

«  rf^  » 


+  0(al) 


(16) 


This  can  be  solved  through  standard  techniques  to  yield  the  time  evolving  variances.  It  must 
be  noted  that  by  solving  only  these  three  equations,  the  order  of  the  solution  has  been  reduced 
to  ): 


where 


/  <<  c* »  ^ 

«C’7» 

\  ^ 

1/Ai 

1/2 

1/Ai 


-I/2A2  -I/A3 

-1/4  -1/4 

-I/A2  -I/2A3 


\  /  \ 

+  0  (17) 


and 


Ai  =  -7 

A2  =  — 7  +  2»t(; 


A2  =  — 7  —  2iw 


If  one  is  interested  in  the  transient  response  and  desires  a  solution  of  order  0{a^ ),  one  could 
cast  all  of  Equations  8  to  15  into  a  matrix  equation  of  the  form  of  Equation  16.  This  would  yield 
a  matrix  of  rank  8.  To  solve  this  analytically  would  require  the  analytical  eigensolution  of  this 
rank  8  ma-triv-  This  is  a  difficult  proposition  at  best  with  unclear  practical  need.  Instead,  for 
each  particular  set  of  parameters,  7  and  a  numerical  matrix  is  obtained.  The  eigensolution 
is  then  obtained  numerically  and  the  results  calculated.  It  should  be  noted  that  there  is  no 
theoretical  loss  of  accuracy  in  solving  the  system  of  equations  this  way;  only  roundoff  error 
degrades  the  accuracy  of  the  solution. 

Thus  we  have  a  method  for  obtaining  in  closed  form  the  time  evolving  moments  of  the 
Duffing  oscillator  subjected  to  a  white  noise  forcing  function. 
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3  Results 

The  response  of  a  Duffing  oscillator  with  damping  coefficient  =  1.0  excited  by  white  noise 
of  intensity  a  =  0.1,  was  calculated  using  Equation  17  and  by  solution  of  all  of  Equations  8 
to  15.  As  a  point  of  comparison,  a  Monte  Carlo  experiment  simulating  a  Duffing  oscillator 
with  the  same  parameters  was  performed.  This  Monte  Carlo  experiment  consisted  of  one 
thousand  iterations  of  a  fourth  order  Runge-Kutta  integration  of  the  following  restatement  of 
Equations  3: 

=  «1  -^V  =  -  V  -  X  -  x^.  (18) 

The  results  of  the  order  0(a3)  analysis,  the  order  0(a2)  analysis,  and  the  Monte  Carlo 
experiment  are  plotted  in  Figures  1  to  3. 

Figure  1  shows  the  time  evolution  of  <  >  as  calculated  by  all  three  methods.  The 

order  0(a3)  analysis  shows  <  >  increasing  monotonically  to  its  steady  state  value  of 

approximately  one.  The  higher  order  analysis  shows  <  >  increatsing  in  nearly  monotonic 

fashion  to  its  steady  state  value  of  about  0.75.  However  this  curve  does  exhibit  some  overshoot 
at  about  2.5  seconds.  The  results  of  the  Monte  Carlo  analysis  are  quite  close  to  those  of  the 
higher  order  analysis,  reaching  the  same  steady  state  value.  The  main  difference  between  the 
two  curves  is  the  slightly  greater  rise  time  of  the  Monte  Carlo  results.  The  higher  limit  of 
the  order  C7(ai)  analysis  and  the  overshoot  of  the  order  C7(a^)  analysis  are  artifacts  of  the 
expansion  process.  The  higher  order,  time  derivatives  differ  from  the  lower  order, 

0(ai),  ones  by  the  subtraction  of  a  <  >  in  the  case  of  ^  <  Cv  >  2a  <  >  in  the 

case  of  ^  <  C*  >•  However  it  can  be  seen  from  Figure  4  that  <  >  >  0  during  the  entire 

time  span  of  interest  and  rises  in  magnitude  quickly,  and  <  >  is  greater  than  zero  during 

the  first  half  of  the  time  span  of  interest  where  the  second  order  moments  are  changing  most 
rapidly.  This  would  explain  the  slower  rise  in  the  magnitude  of  the  order  C7(at )  solutions  than 
the  order  0(a3  )  solutions.  It  also  indicates  that  inclusion  of  the  higher  order  terms  serve  to 
lower  the  overall  value  of  the  analytical  results  and  that  the  omission  thereof  causes  the  higher 
overall  values  of  the  analytical  methods  than  those  of  the  Monte-Carlo  methods.  As  one  would 
expect,  the  higher  order  analyses  consistently  show  mpre  points  of  inflection  than  do  the  lower 
order  analyses. 

Figure  2  shows  nearly  identical  curves  for  all  methods.  As  in  Figure  1,  the  order  C7(at) 
analysis  shows  more  overshoot  in  one  place  than  does  the  lower  order  analysis.  Here  this 
overshoot  occurs  at  about  four  seconds,  at  the  second  local  extreme.  It  is  also  noted  that, 
only  the  order  C7(al)  method  accurately  reflects  the  re^on  where  the  Monte-Carlo  curve  is 
negative.  The  initial  excursion  of  the  Monte  Carlo  curve  is  not  as  great  as  that  of  the  other 
two  curves,  although  all  three  approach  the  x-axis  as  time  increases. 

Figure  3  Provides  a  qualitatively  similar  comparison  of  results  with  both  analyses  almost 
coincidental  and  again  slightly  greater  in  the  transient  region  than  the  Monte  Carlo  experiment. 
All  three  curves  approach  a  steady  state  value  of  approximately  1. 

Figure  5  shows  how  the  time  evolution  of  <  >,  as  calculated  by  the  order  0(a3 )  analysis, 

is  affected  by  increasing  7.  Figure  6  shows  the  same  study  for  the  Monte  Carlo  experiment. 
Figure  5  shows  not  only  a  decrease  in  the  steady  state  values  of  <  (^  >  with  increasing  values 
of  7,  but  also  a  smoothing  of  the  curves.  The  decrease  in  the  steady  state  values  is  also  seen 
in  the  Monte  Carlo  curves.  The  cause  of  the  decrease  is  a  physical  one:  increased  damping 
implies  increased  viscous  energy  dissipation  which  leads  to  smaller  excursions  of  the  oscillator. 
Comparison  of  Figures  5  and  6  shows  that  even  when  7  becomes  large  enough  to  violate 
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the  order  unity  assumption,  the  analysis  still  gives  reasonable  results.  In  fact  the  agreement 
between  results  becomes  better  at  higher  values  of  7.  It  can  be  seen  that,  even  for  7  =  0.6, 
the  Monte  Carlo  curves  are  relatively  smooth  with  the  exception  of  the  small  scale  “wiggles” 
inherent  in  Monte-Carlo  analysis.  This  indicates  that  the  multiple  local  maxima  and  minima, 
i.e.,  at  3,  5,  7,  9  seconds,  in  the  analytical  curves  at  7  =  0.6  are  due  to  the  low  order  of  7  of  the 
analysis.  It  was  assumed  at  the  beginning  of  the  analysis  that  7  is  of  order  unity.  However  it 
can  be  seen  from  these  curves  that  the  further  7  is  from  this  assumption,  the  worse  the  results 
of  the  analysis. 

The  behavior  of  the  system  as  7  becomes  small  is  further  investigated  in  Figures  7  and  8. 
These  curves  show  how,  as  7  goes  from  0.4  to  0.3  the  analysis  breaks  down  completely.  It 
was  shown  in  Figures  5  and  6  that  the  7  =  0.3  analysis  differed  signiiicamtly  from  the  Monte 
Carlo  analysis  in  the  transient  region  of  0  to  5  seconds,  although  its  large  time  behaviour 
was  accurate.  Figures  7  and  8  show  that  as  7  is  made  smaller,  even  the  large  time  behaviour 
becomes  unreasonable  with  <  >  becoming  negative  which  is  clearly  impossible. 

Figures  9  and  10  show  how  the  time  evolutions  of  '<  >  are  affected  by  increasing  values 
of  a.  It  can  be  seen  in  Figure  9  what  happens  as  the  assumption  that  a  is  small  as  compared  to 
unity  is  violated.  As  is  shown  in  Figure  9,  the  curve  for  a  =  0.3  exhibits  oscillations  that  damp 
out  slowly  with  time.  As  a  becomes  even  larger,  negative  values  for  <  >  develop.  However, 

since  the  average  of  a  squared  quantity  cannot  be  negative,  these  results  are  spurious.  The 
trend  of  the  steady  state  value  of  <  >  decreasing  with  increasing  values  of  a  is  common 

to  both  figures.  This  does  not  violate  the  basic  assumption  that  (  is  of  order  unity  when  a 
is  small.  (  does  indeed  remain  of  order  unity;  it  was  not  assumed  that  C  was  independent 
of  a.  <  >  decreases  with  increasing  a  due  to  the  nonlinearity  of  the  Duffing  oscillator. 

The  energy  stored  in  the  nonlinear  spring  is  greater  than  that  stored  in  the  linear  spring  by 
^noniinear  -  ^(tnear  =  Therefore  the  effect  of  the  nonlinearity  of  the  Duffing  oscillator  is  to 
decrease  <  >,  and  this  effect  will  increase  with  increasing  <  >  and  therefore  increasing 

<  a  > 

Figures  11  to  13  show  the  results  of  Monte-Carlo  experiments  for  the  system  driven  by 
colored  noise.  In  each  curve  the  time  evolving  behaviour  of  the  second  order  moment  is 
depicted  for  various  values  of  the  correlation  time,  Vg.  All  three  figures  show  essentially  the 
same  behaviour  for  the  case  of  Tg  =  0.001  and  rg  s  0.01  as  for  the  white  noise  case.  Figures  1 
to  3.  For  the  case  of  Tg  =  0.1,  the  of  the  correlated  nature  of  the  noise  is  noticeable,  but 
perhaps  acceptable  for  some  applications.  It  is  clear  that  when  the  correlation  time  becomes 
greater  than  0.1,  the  results  differ  significantly  from  the  white  noise  case.  This  is  as  one  would 
expect:  at  Tg  >  0.1  the  correlation  time  becomes  comparable  to  the  natural  period  of  the 
oscillator,  which  is  about  one  second.  The  correlated  nature  of  the  noise  appears  as  a  effect 
of  time  scale  Tg  on  the  time  history  of  the  correlated  noise.  If  the  time  scale  of  the  correlation 
is  much  smaller  than  the  natural  period  of  the  osdllator,  then  the  oscillator  caimot  respond 
to  this  effect.  However  as  Tg  approaches  the  natural  period  of  the  oscillator,  the  oscillator  can 
be,  and  is,  affected. 
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Figure  4:  The  time  evolution  of  the  fourth  order  moments  for  7  =  1,0,  a  =  0.1. 
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Figtire  5:  <  >  vs  Time  for  increasing  values  of  7,  a  =  0.1,  calculated  analytically  to  order  0{a\) 
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Figure  8:  <  (*  >  versus  T  for  decreasing  values  of  7,  o  =  0.1,  calculated  by  Monte-Carlo  simulation 
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Figure  11: 
simulation 


Figure  12: 
simulation 


<  C*  >  versus  T  for  different  values  of  Tc,  of  7  s=  1.0,  a  =  .1,  Calculated  by  Monte-Carlo 


<  Cv  >  versus  T  for  different  values  of  t^,  of  7  =  1.0,  a  =  .1,  Calculated  by  Monte-Carlo 


Time,  seconds 


Figure  13:  <  rj^  >  versus  T  for  different  values  of  Tc,  7  =  1.0,  a  =  .1,  Calculated  by  Monte-Carlo 
simulation 


4  Conclusions 


It  can  be  seen  from  the  figures  that  the  agreement  between  the  results  of  the  various  methods 
of  solution  is  good.  The  trends  observed  are  not  surprising.  The  further  the  parameters  of 
the  system  are  from  violating  the  assumptions  of  the  ^alysis,  the  better  the  analytical  results 
agreed  with  those  of  the  Monte  Carlo  experiment.  It  is  also  seen  that  by  and  large,  the 
higher  order  analysis  is  more  accurate.  One  interesting  point  is  that  the  analytical  techniques 
are  consistently,  albeit  slightly,  greater  than  the  Monte  Carlo  experiments.  This  conservative 
nature  of  the  analytical  method  should  be  noted  in  applying  it. 

The  method  can  be  applied  to  systems  where  the  three  basic  conditions  of  the  expansion  are 
met:  that  the  intensity  of  the  forcing  fimction,  a  be  of  order  smaller  than  unity,  the  damping 
of  the  system  be  of  order  at  least  unity,  and  the  forcing  function  be  essentially  uncorrelated  in 
time,  i.e.,  correlation  time  <  0.1 

This  technique  is  a  vafid  analytical  tool.  It  is  well  suited  to  studying  the  behavior  of  a 
system  under  a  variety  of  conditions.  The  results  it  gives  are  accurate  and  computationally 
fast  enough  to  embed  such  a  model  as  an  element  of  a  larger  model. 
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Abstract 


In  Rodriguez  and  Van  Kampen’s  1976  paper  [1],  a  method  of  extracting  information  from 
the  Fokker-Planck  equation  without  having  to  solve  the  equation  is  outlined.  The  Fokker- 
Planck  equation  for  a  Duffing  oscillator  excited  by  white  noise  is  expanded  about  the  intensity 
of  the  forcing  function,  a.  In  Wdnstein  and  Benaroya,  the  effect  of  the  order  of  expansion  is 
investigated  by  carrying  the  expansion  to  a  higher  order.  The  effects  of  varying  the  system 
parsoneters  is  also  investigated.  All  results  are  verified  by  comparison  to  Monte  Carlo  exper¬ 
iments.  In  this  paper,  the  van  Kampen  expansion  is  modified  and  applied  to  the  case  of  a 
Duffing  osdllator  excited  by  colored  noise.  The  effect  of  the  correlation  time  is  investigated. 
Again  the  results  are  compared  to  those  of  Monte  Carlo  experiments.  It  was  found  that  the 
expansion  compared  closely  with  those  of  the  Monte  Carlo  experiments  as  the  correlation  time, 
Te,  was  varied  from  0.001  to  10  seconds.  Examination  of  the  results  revealed  that  the  colored 
noise  can  be  categorized  in  one  of  four  ways:  1)  for  Tc  <  Cl(0.01)sec  the  noise  can  be  considered 
as  white  for  all  intents  and  purposes,  2)  for  =  0(0.1)  the  noise  could  be  considered  white 
for  some  purposes,  3)  for  Tg  =  O(l.Osec)  the  correlated  nature  of  the  noise  must  be  considered 
in  an  analysis,  and  4)  for  0(1.0)  <  Tg  the  noise  can  be  considered  as  deterministic. 
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1  Introduction 


The  Fokker-Planck  equation  has  proven  to  be  a  useful  tool  in  the  analysis  of  simple  nonlinear 
osdllators  excited  by  stochastic  processes.  As  a  partial  differential  equation  for  the  probability 
density  function  of  the  response,  its  solution  completely  defines  the  solution  of  the  problem. 
It  can  be  used  to  analyze  both  a  single  oscillator  of  the  form 

mx +  'i{x,x)x  +  k(x,x)x  =  T{t)y  (1) 

or  a  system  of  multiple,  linked,  oscillators  of  the  form 

Mi  +  r(i,  x)x  +  K(i.  x)x  =  £(t).  (2) 

In  many  cases,  a  physical  system  can  be  approximated  by  such  a  system  of  nonlinear 
oscillators.  The  systems  so  modeled  can  range  hrom  a  Brownian  particle  to  structures  excited 
by  von  Karmann  vortex  shedding.  Such  modeUng  can  be  useful  for  gaining  insight  into  a 
problem  and  the  way  in  which  the  system  will  behave  as  certain  parameters  are  varied. 

Once  one  has  decided  on  the  system  of  oscillators  to  be  used  to  represent  the  physical 
system,  the  derivation  of  the  Fokker-Planck  equation  is  relatively  straightforward,  although 
tedious.  The  problem  remains  of  how  to  solve  it  for  the  probability  distribution  of  the  response. 
In  a  very  few  cases,  the  Fokker-Planck  equation  can  be  solved  analytically,  but  in  most  cases 
no  analytical  solution  exists  and  one  usually  must  resort  to  a  niunerical  solution.  However  this 
can  be  computationally  intensive  and  ^ves  little  insight  into  the  larger  problem. 

In  their  1976  paper,  Eodiiguez  and  van  Kampen  outline  a  method  of  dealing  with  the  case 
of  an  oscillator  excited  by  weak  Gaussian  white  noise.  The  Fokker-Planck  equation  of  the 
system  is  expanded  about  the  intensity,  a,  of  the  driving  function.  This  expansion  is  carried 
to  the  order  O(a^).  In  this  way  the  statistics  of  the  fluctuations  are  obtained  directly.  This 
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I 


I 


I 
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method  shows  promise  as  a  way  to  use  the  Fokker-Planck  equation  to  gain  useful  information 
about  a  v/idet  variety  of  systems  than  was  possible  before. 

This  is  the  second  of  a  planned  series  of  papers  exploring  the  usefulness  of  this  method.  In 
the  previous  paper  [2],  the  method  was  applied  to  the  problem  of  a  Duffing  oscillator  excited 
by  Gaussian  white  noise.  The  inherent  assumptions  of  the  method  were  explained  there  in 
detail,  the  expansion  was  carried  both  to  the  same  order  as  in  the  original  paper,  and  to  order 
C7(al),  and  results  were  presented  and  compared  to  Monte  Carlo  experiments.  Parametric 
studies  were  also  performed  on  the  parameter  of  expansion  as  well  as  on  the  other  important 
variable  in  the  expansion:  the  coefficient  of  damping. 

In  this  paper,  the  expansion  is  applied  to  the  case  of  a  Duffing  oscillator  excited  by  expo¬ 
nentially  correlated  noise.  A  parametric  study  is  performed  on  the  correlation  time,  and 
the  results  are  compared  to  those  of  Monte  Carlo  experiments.  The  range  of  Tc  for  which  the 
correlated  noise  can  be  treated  as  white  is  identified  as  well  as  the  range  of  Tc  for  which  the 
correlated  nature  of  the  noise  cannot  be  ignored. 
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2  Expansion  of  the  Fokker  Planck  Equation  for  a 
Duffing  Oscillator  Excited  By  Correlated  Noise 

In  Gang  [3],  a  method  is  outlined  for  deriving  the  Fokker  Planck  Equation  of  a  system  driven 
by  colored  noise.  The  procedure  outlined  consists  of  defining  a  dummy  variable,  y,  to  represent 
the  colored  noise  process,  y  is  defined  as  a  first  order,  linear,  dififerential  filter  of  the  white 
noise  process: 

ciy{t)  =  C2y(t)  +  czTit),  (3) 

where  Ci  are  constants  and  is  the  Gaussian  white  noise  process  defined  by 

<  F{t)  >  =  0 

<F(t)F(t')>  =  2aS(t-t').  (4) 

The  system  is  now  formulated  as  a  system  of  first  order  differential  equations  driven  by 

«(1): 

i  =  (5) 

Cl  Cl 

X  =  F*(x,t>,y)  (6) 

V  =  F„(x,w,y).  (7) 

Several  examples  of  colored  noise  filters  enst  in  the  literature.  Billah  and  Shinozuka  [4] 
use  the  following 

rcy(t)  =  -y(t)  +  F(t),  (8) 

where  Tc  is  the  correlation  time.  It  can  be  seen  that  as  Tc  0,  y(t)  F(t).  For  this  derivation, 
it  will  be  assumed  that  Tc  is  not  vanishingly  small,  so  that  Equation  8  can  be  written  as: 

y(t)  =  -fiy(t)  +  fiF(t),  (9) 
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where  /*  = 

The  Duffing  oscillator  excited  by  colored  noise  can  now  be  written  as 

*(0  +  7*(0  +  *(0  +  =  »(0.  (10) 

where  7  is  the  coefficient  of  damping. 

By  defining  v(t)  =  i(t),  the  system  can  be  recast  as  three  linked,  ordinary  differential 
equations  in  time  in  the  form  of  Equations  5  to  7: 


x(t)  = 

v(t) 

(11) 

v(t)  = 

y(t)  -  7t;(t)  -  z(t)  -  x®(t) 

(12) 

y(t)  = 

-M*)  +  M^(0* 

(13) 

Define  /(z,  u,  y\  t)  to  be  the  joint  probability  density  function  of  i(t),  t;(t),  and  y{t)  at  time 
t.  The  Fokker- Planck  equation  for  f(x,v,y\t),  can  now  be  derived  as 

Equation  14  is  the  governing  equation  for  the  time  evolution  of  the  transition  probability 
density  function  f(x,  v,  y;  t).  EYom  this  point  the  derivation  of  the  previous  paper  [2]  will  be 
followed. 

As  was  shown  in  the  previous  paper,  the  response  of  the  oscillator  can  be  separated  into  a 
deterministic  component  due  to  the  initial  conditions  and  a  random  component  of  magnitude 
However,  by  assuming  that  the  oscillator  is  initially  at  rest,  the  deterministic  com¬ 
ponent  can  be  shown  to  be  equal  to  zero.  Therefore  the  following  substitutions  are  made  into 
Equation  14. 

X  =  y/aC  (15) 
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t;  =  y/arj 

(16) 

II 

(17) 

f  {y/aC,\^n,Vop\t)  =  a  ®/^n(C,  »?,/>;<)• 

(18) 

The  factor  will  be  omitted  from  the  defimtion  of  n(C,T;,p;t).  If  carried  through  the 

derivations,  it  would  be  divided  out  at  a  later  stage  of  the  derivation. 

It  must  be  explained  why,  as  is  implied  in  Equation  17,  y  is  of  the  order  of  magnitude  . 
This  is  because  the  magnitude  y  is  of  the  same  order  of  magnitude  of  the  white  noise  forcing 
function,  F.  But  it  can  be  seen  from  the  definition  of  f',  Eqs  4,  that  a  =  2ap.  Therefore  F, 
and  consequently  y,  are  of  0(q^  ). 

The  relationships  between  the  partial  derivatives  of  /  and  n  are  obtained  as: 


a 

a 


an 

dc 

an 

dv 


idf  an 

“’ay  "  dp 

^  _  an 
at  at 


(19) 

(20) 
(21) 
(22) 


This  yields  the  following  transformed  FPK: 

By  manipulating  the  left  hand  side  of  Equattion  23,  the  following  is  obtained: 


an.  an.  an. 
ac^  a»7^  ap^ 


+^^(pn)+,.^^n. 


(24) 


If  one  substitutes  the  definitions  of  *7»  ^<1  P  into  Eqs  11,  12  and  13,  one  can  use  the 


resultant  equations  to  separate  Equation  24  into  the  following  three  equations: 


cn  =  i7n 


(25) 
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(26) 


j)n  =  {2p  -  7»?  -  c  -  aC®}  n. 

pVL  =  -ppH  -  (27) 

Using  Eqs  27,  25,  and  26,  one  can  find  the  time  derivatives  of  the  higher  order  moments  of 
and  p.  As  an  example,  the  time  derivative  of  <  >  ^jji  derived.  To  derive  ^  <  P^  >■, 

one  multiplies  Equation  27  by  2p,  notes  that  2pp  =  and  gets 

p^n  =  -pp^n  -  p^pllp.  (28) 

Integrating  this  over  the  range  of  all  three  arguments  of  n  and  noting  the  definition  of  expec¬ 


tation,  it  is  found: 

^  <  p^  >=  -2p  <  p^  >  -2p2  J  j  J ^pUp  dp  dC  dr].  (29) 

The  second  term  on  the  right  hand  side  can  be  expanded  as  follows: 

/ -  j  j  j  !('’•')»  -O]  dplt(  (fi, 

= -i+//”i/>n);:r„  =  -i.  (30) 

This  last  step  is  possible  because  the  existence  of  <  p  >  guarantees  that  pll  — »  0  as  p  -»  ±oo. 
Equation  29  becomes: 

^  <  p2  >=  _2p  <  p2  >  -|-2p^  (31) 

The  time  derivatives  of  the  other  second  order  moments  are  similarly  obtained: 

^<C'>  =  2<«>=-2<Ci7>  (32) 

=  2<'rt>=-2<C0> 

-27  <  >  -1-4  <  qp  >  -2a  <  >  (33) 

^  <  C»?  >  =  <  C*?  >  +  <  >=  -  <  >  -7  <  C»7  > 
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• 

-  <  17^  >  +2  <  Cp  >  -a  <  C*  > 

(34) 

2 

-2p  <  p^  >  -2p^ 

(35) 

• 

II 

A 

V 

-H  <Cp>  -  <rip> 

(36) 

^<v>  = 

-  <  Cp  >  -7  <  i/p  >  +2  <  p*  >  -o  <  C^p  >  . 

(37) 
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3  Results  and  figures 


Figures  1,  3  and  5  show  the  results  of  the  expansion  for  a  =  0.1,  7  =  1.0,  and  different  values 
of  the  correlation  time,  Tc'.  Tc  =  0.001,  Tc  =  0.01,  Tc  =  0.1,  Tc  =  1.0,  Tc  =  10.0.  Figures  2,  4 
and  6  show  the  results  of  Monte  Carlo  experiments  for  similar  cases. 

As  in  the  previous  paper,  there  is  very  good  agreement  between  the  two  methods  with 
similar  differences  as  well.  The  anal}^tical  methods  consistently  give  results  slightly  greater 
magnitude  than  do  the  Monte  Carlo  experiments.  This  phenomenon  was  seen  in  the  previous 
paper.  The  anal)rtical  results  are  also  consistently  smoother,  showing  none  of  the  small  time 
scale  fluctuations  inherent  in  the  Monte  Carlo  technique. 

It  can  be  seen  in  all  the  ctirves  that  the  traces  representing  =  0.001  and  Tc  =  0.01  are 
almost  identical.  In  some  cases,  most  notably  Figures  1  and  3,  the  two  traces  are  almost 
indistinguishable.  The  difference  between  the  traces  representing  Tc  =  0.001  and  Tc  =  0.01 
are  slightly  more  pronounced  in  the  Monte  Carlo  results.  However,  the  difference  between  the 
Tc  =  0.001  and  Tc  =  0.01  traces  of  the  Monte  Carlo  results  is  only  of  the  order  of  the  small 
time  scale  fluctuations  inherent  in  the  Monte  Carlo  technique.  These  figures  indicate  that  for 
Tc  of  order  0(0.01)  or  less,  the  noise  can  be  assumed  to  be  uncorrelated,  or  white.  The  traces 
representing  the  response  for  rc  =  0.1  differ  noticeably  from  those  representing  the  results  for 
Tc  <  0.1.  However,  even  for  Tc  =  0.1,  the  results  are  still  quite  close  to  those  for  Tc  <  0.1  and 
the  white  noise  approximation  may  stUl  be  useful  for  some  uses.  One  would  expect  the  effects 
of  the  correlated  nature  of  the  noise  to  become  significant  at  about  Tc  =  0.1;  at  Tc  =  0.1,  the 
correlation  time  begins  to  become  comparable  to  the  natural  period  of  the  oscillator,  which 
is  about  one  second.  The  correlated  nature  of  the  noise  appears  as  an  effect  of  time  scale, 
Tc  on  the  time  history  of  the  correlated  noise.  If  the  time  scale  of  the  correlation  is  much 


smaller  than  the  natural  period  of  the  oscillator,  then  the  oscillator  cannot  respond  to  this 
effect.  However,  as  Tc  approaches  the  natural  period  of  the  oscillator,  the  oscillator  can  be, 
and  is,  affected. 

It  can  likewise  be  seen  that  as  the  correlation  time  of  the  noise  becomes  much  greater 
than  the  natural  period  of  the  oscillator,  the  m^nitude  of  the  random  response  approaches 
zero.  This  is  because  when  the  correlation  time  is  much  larger  than  the  natural  period  of  the 
oscillator,  the  oscillator  responds 
show  little  random  response  for 


to  the  noise  as  if  it  were  deterministic.  Hence  the  figures 
=  1.0  seconds  and  almost  none  at  all  for  Tc  =  10.0  seconds. 
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Figure  1:  <  >  vs  Time  for  different  values  of  7  =  1.0,  a  =  0.1,  calculated  analytically  to 

order  0{a^). 


Unie,  seconds 


Figure  2:  <  C*  >  vs  Time  for  different  values  of  Tc,  7  =  1.0,  a  =  0.1,  calculated  by  Monte-Carlo 


simulation. 


Figure  3:  <  C»?  > 
order  0{a^). 


Figure  4:  <  Cv  > 
simulation. 
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vs  Time  for  different  values  of  7  =  1.0,  a  =  0.1,  calculated  analytically  to 


Time,  seconds 


vs  Time  for  different  values  of  Tg,  7  =  1.0,  a  =  0.1,  calculated  by  Monte-Carlo 
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Figure  5:  <  r)^  >  vs  Time  for  different  values  of  Tc,  7  =  1.0,  a  =  0.1,  calculated  analytically  to 
order  0(at). 


Time,  seconds 


Figure  6:  <  >  vs  Time  for  different  values  of  Tc,  7  =  1.0,  a  =  0.1,  calculated  by  Monte-Carlo 


simulation. 


4  Conclusions 


The  overwhelming  similarity  between  the  results  ^ven  by  the  two  methods  implies  that  this 
adaptation  of  the  van  Kampen  expansion  is  an  accurate  tool  for  predicting  the  statistics  of  the 
response  of  an  oscillator  excited  by  colored  noise.  However,  it  was  also  seen  that,  depending 
on  the  magnitude  of  the  correlation  time,  perhaps  as  compared  to  the  natural  period  of  the 
oscillator,  simplifying  assumtions  can  be  made  that  obviate  the  need  for  this  adaptation. 
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The  van  Kampen  expansion  for  a  linked  Duifing-linear 
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Abstract 

In  Rodrignez  and  Van  Kampen ’s  1976  paper  [1],  a  method  of  extracting  information  from 
the  Fokker-Planck  equation  without  having  to  solve  the  equation  is  outlined.  The  Fokker- 
Planck  equation  for  a  Duffing  oscillator  excited  by  white  noise  is  expanded  about  the  intensity 
of  the  forcing  function,  a.  In  Weinstdn  and  Benaroya  [2],  the  effect  of  the  order  of  expansion 
is  investigated  by  carrying  the  expansion  to  a  h^her  order.  The  effects  of  varying  the  system 
parameters  is  also  investigated.  All  results  are  verified  by  comparison  to  Monte  Carlo  experi¬ 
ments.  In  Weinstein  and  Benaroya  [3]  the  van  Kampen  expansion  is  modified  and  applied  to 
the  case  of  a  Duffing  oscillator  excited  by  colored  noise.  The  effect  of  the  correlation  time  is 
investigated.  Again  the  results  are  compared  to  those  of  Monte  Carlo  experiments.  In  both 
cases,  the  results  of  the  analyses  agreed  closely  with  those  of  the  Monte  Carlo  experiments. 
In  this  paper  the  van  Kampen  expansion  is  applied  to  linked  linear-Duffing  oscillators.  Again, 
parametric  studies  are  done  on  the  system  parameters  and  the  correlation  coefficient  of  the 
driving  force  and  the  results  compared  to  those  of  Monte  Carlo  experiments.  It  was  found 
that  the  analytical  results  compared  closely  with  the  numerical  results  as  long  as  the  initial 
assumptions  of  the  expansion  are  not  violated. 


1  Introduction 

The  Fokker-Planck  equation  has  proven  to  be  a  useful  tool  in  the  analysis  of  simple  nonlinear 
oscillators  excited  by  stochastic  processes.  As  a  partial  differential  equation  for  the  probability 
density  function  of  the  response,  its  solution  completely  defines  the  solution  of  the  problem. 

It  can  be  used  to  analyze  both  a  single  oscillator  of  the  form 

mx  +  7(i,  x)x  +  k(x,  x)x  =  ^(t),  (1) 

or  a  system  of  multiple,  linked,  <M3cillators  of  the  form 

Mx  +  r(x,  x)i  -I-  K(i,  x)x  =  Z(t).  (2) 
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In  many  cases,  a  physical  system  can  be  approximated  by  such  a  system  of  nonlinear 
oscillators.  The  systems  so  modeled  can  range  from  a  Brownian  particle  to  structures  excited 
by  von  Karmann  vortex  shedding.  Such  modeling  can  be  useful  for  gaining  insight  into  a 
problem  and  the  way  in  which  the  system  will  behave  as  certain  parameters  are  varied. 

Once  one  has  decided  on  the  system  of  oscillators  to  be  used  to  represent  the  physical 
system,  the  derivation  of  the  Fokker-Planck  equation  is  relatively  straightforward,  although 
tedious.  The  problem  of  how  to  solve  it  for  the  probability  distribution  of  the  response  remains. 
In  a  very  few  cases,  the  Fokker-Planck  equation  can  be  solved  analytically,  but  in  most  cases 
no  analytical  solution  exists  and  one  usually  must  resort  to  a  numerical  solution.  However, 
this  can  be  computationally  intensive  and  ^ves  little  insight  into  the  larger  problem. 

In  their  1976  paper,  Rodriguez  and  van  Kampen  outline  a  method  of  dealing  with  the  case 
of  an  oscillator  excited  by  weak  Gaussian  white  noise.  The  Fokker-Planck  equation  of  the 
system  is  expanded  about  the  intensity,  o,  of  the  driving  function.  This  expansion  is  carried 
to  the  order  C7(a^).  In  this  way  the  statistics  of  the  fluctuations  are  obtained  directly.  This 
method  shows  promise  as  a  way  to  use  the  Fokker-Planck  equation  to  gmn  useful  information 
about  a  wider  variety  of  systems  than  was  possible  before. 

This  is  the  of  a  series  of  papers  exploring  the  usefulness  of  this  method.  In  this  paper, 
the  van  Kampen  expansion  is  applied  to  coupled  linear-Duffing  oscillators.  This  is  done  for 
two  reasons.  The  first  is  to  demonstrate  that  the  expansion  can  be  applied  to  a  system 
of  coupled  osdllators.  The  second  is  that  some  physical  systems  are  well  modeled  by  such 
coupled  oscillators.  An  example  of  such  a  system  is  an  offshore  structure  under  the  influence 
of  ocean  waves  [4,  5].  These  waves  have  a  certain  periodic  nature,  but  are  of  random  height 
and  intensity.  The  structure  itself  can  be  modelled  as  a  series  of  oscillators. 

There  are  several  ways  in  which  such  coupled  oscillators  can  be  formulated.  The  external 
force  can  drive  either  the  linear  or  the  Dufling  oscillator.  Here,  it  is  decided  to  have  the 
external  force  drive  the  linear  oscillator.  The  formulation  of  the  solution  would  be  the  same 
if  the  situation  was  reversed.  There  are  also  several  ways  in  which  the  two  oscillators  can  be 
coupled:  through  a  damping  function,  through  a  spring  function,  or  through  a  combination  of 
the  two.  In  this  analysis,  the  coupling  chosen  is  of  the  form  of  a  simple  linear  spring.  Here, 
the  system  is  formulated  in  such  a  way  that  the  effect  of  the  Duffing  oscillator  would  not 
feed  back  into  the  linear  oscillator.  Such  a  system  would  more  accurately  model  a  system 
where  the  driving  force  is  an  overwhelming  physical  phenomenon  impervious  to  the  effect  of 
the  structure  being  modeled.  However,  it  appears  that  all  of  the  above  cases  can  be  modeled 
using  the  following  technique. 


2  Derivation  and  Expansion  of  the  Fokker-Planck 
equation 


The  first  step  in  the  derivation  of  the  Fokker-Planck  equation  for  any  system  is  the  formu¬ 
lation  of  the  governing  equations.  The  governing  equations  for  the  system  described  in  the 
Introduction  are: 


2/(<)  +  7^(0 
£i(t)  +  7iii(0  +  ®i(0 
xa(t)  +  73®2(f)  +  *2(0  +  £*3(0 


(3) 

y(<) 

(4) 

k  [3;i(t)  -  X2(t)] , 

(5) 

2 


where  F{t)  is  defined  by 


<Fit)>  =  0, 

<F{t)F{i^)>  =  2a6{t-t').  (6) 

That  the  solution  of  Ekjuation  3  is  exponentially  colored  noise  is  shown  in  several  sources  such 
as  Billah  and  Shinozuka  [6]. 

Because  the  Fokker-Planck  equation  requires  that  the  governing  equations  be  cast  as  a 
series  of  first  order  difierential  equations,  the  following  new  variables  are  defined: 


vi  =  ii  (7) 

V3  =  *2.  (8) 


Using  these  variables  the  governing  equations  can  be  rewritten  in  the  following,  equivalent, 
form: 


y(t)  = 

(9) 

ii(t)  = 

Vi(t) 

(10) 

Vx{t)  = 

y(t)  -  iiviit)  -  xi(t) 

(11) 

*2(0  = 

V2(t) 

(12) 

V2it)  = 

kxi{t)  -  72«2(t)  -  (1  +  k)X2it)  -  £xf(t). 

(13) 

The  Fokker-Planck  equation  of  this  system  can  be  derived  as: 


-  72»3  -  (1  +  k)x2  -  €xl)f]  + 


(14) 


As  was  shown  in  Weinstein  and  Benaroya  [2],  the  response  of  the  oscillator  can  be  separated 
into  a  deterministic  component  arising  from  the  initial  conditions,  and  a  random  component  of 
magnitude  0{y/^.  However,  by  assuming  the  osdllator  to  be  initially  at  rest,  the  deterministic 
component  can  be  shown  to  be  equal  to  zero.  Therefore,  the  following  substitutions  are  made 
into  Equation  14: 


y  =  y/ap 
®i  =  VaCi 

vi  =  y/a  rti 
X2  =  Cz 
V2  =  y/arfy 

f{y,xi,vi,x2,v2;t)  =  Q!~fn(p,Ci,»7i.C2,»?2;t). 


(15) 

(16) 

(17) 

(18) 

(19) 

(20) 


n(p,Cii’7i)C2>*72if)  is  til®  joint  probability  density  of  the  transformed  variables.  The  factor 
a~l  will  be  omitted  from  the  definition  of  H.  If  carried  through  the  derivations,  it  would  be 
divided  out  at  a  later  stage. 
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The  relationships  between  the  partial  derivatives  of  /  and  n  are  easily  obtained  as: 

(21) 
(22) 

(23) 

(24) 

(25) 

(26) 

Equations  15  to  19,  as  well  as  the  equations  above,  Eqs  21  to  26,  can  be  substituted  into 
the  Fokker-Planck  equation.  Equation  14,  yielding: 

Ht  =  /I  (pU\  -  -  [(p  -  7iih  -  Ci)n],n  - 

-  [(*Ci  -  72»?j  -  (1  +  fc)C2  -  (27) 

The  left  hand  side  of  Equation  27  is  transformed  into  partial  derivatives  of  11  in  the  trans¬ 
formed  variables,  p,  Cii  Cii  72i  yielding 

pDp  -  CiDci  -  = 

M  (pD)^  -  fh^Ci  -  l(p  -  7i»h  -  Ci)n],;,  -  %n(j 

-  [(Ki  -  72%  -  (1  +  *)C2  -  ear]^)Tl]^  + 

(28) 

One  can  now  substitute  the  definitions  of  P,  Cl>  *71  C2>  ^<1  %  Equations  10  to  13  and 
multiply  by  a~^.  The  resulting  four  equations  are: 


Cl 

=  m 

(29) 

m 

=  p-ifim-Ci 

(30) 

Cl 

=  % 

(31) 

m 

=  *Ci -72%- (1 +  • 

(32) 

If  one  multiplies  each  of  the  above  equations  by  11  and  differentiates  Equation  29  with 
respect  to  Equation  30  with  respect  to  »/i,  Equation  31  with  respect  to  C2>  ami  Equation  32 
with  respect  to  ijj,  each  of  the  resulting  equations  can  be  added  to  the  transformed  Fokker- 
Planck  equation.  Equation  28.  The  resultant  is  the  following  equation: 

(33) 

This  can  be  integrated  with  respect  to  p  to  give 

-pn  =  ppn-|-p%.  (34) 
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^  df  an 
^Ty  =  Tp 
r-  df  an 

^  dxx  acx 

r-  df  an 

a«,  drn 

r-  df  ^  ^ 

a*2  a^a 

^  M 

^  dv2  drti 

a/  ^ 

dt  dt " 


1 

► 

Thus,  the  Fokker- Planck  equation  of  the  system  of  linked  oscillators  has  been  transformed 

to  the  following  five  equations: 

pH  =  -ppll  + 

(35) 

CiH  =  i^n 

(36) 

> 

ihn  =  (p-7i»7i -Ci)n 

(37) 

CsD  =  %n 

(38) 

qsll  =  [kCi  -  72»?2  -  (1  +  *)C2  -  eaCi] 

(39) 

As  shown  in  Weinstein  and  Benaroya  [2]^  the  time  derivatives  of  the  second  order  moments  I 

> 

can  be  found  as 

A 

V 

=  2  <  Ci»?i  > 

(40) 

d  2 

=  -271  <  q?  >  -2  <  Cim  >  +2  <  ihp  > 

(41) 

d  ^ 

=  -  <  Cl  >  +  <  »?i  >  -7i  <  Ci»h  >  +  <  CiP  > 

(42) 

=  2  <  C2»?3  > 

(43) 

A 

V 

=  <  »?iC2  >  +  <  C1O2  > 

(44) 

$ 

d 

Jt  <Vl<2> 

=  -  <  C1C2  >  -71  <  V1C2  >  +  <  >  +  <  C2P  > 

(45) 

d  3 

dt<^> 

=  2fc  <  Ci»j2  >  -2(1  +  k)<  C2»?2  >  -272  <  »?2  > 

(46) 

d  ^ 

=  fc  <  Cl  >  -(1  +  fc)  <  C1C2  >  -72  <  Cl»l2  >  +  <  »n*?2  > 

(47) 

=  k  <  Ci»li  >  —(1  +  k)  <  »7iC3  >  —  <  Cl 72  > 

-(71  +  72)  <  mv2  >  +  <  72P  > 

(48) 

=  -(1  +  fc)  <  cl  >  +*  <  C1C2  >  +  <  72  >  -72  <  C272  > 

(49) 

A 

V 

=  -2/i  <  >  -2p? 

(50) 

d  , 

^  <(!(>> 

=  -p  <  Cip  >  +  <  Tip  > 

(51) 

► 

i<mp> 

=  <  P*  >  -  <  Cip  >  -(71  +  p)  <  Tip  > 

(52) 

d  , 
5j<6p> 

=  -p  <  C2P  >  +  <  72P  > 

(53) 

d 

=  fc  <  CiP  >  -(1  +  *)  <  C2P  >  -(72  +  p)  <  72P  >  • 

(54) 

) 

The  previous  15  equations,  Equations  40  to  54,  can  be  cast  as  a  single  matrix  equation. 

=  A2+b, 

(55) 

> 

5 

where 


6 


0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

2 

0 

0 

0 

0 

0 

1 

0 

0 

0 

0 

2 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

1 

0 

0 

0 

0 

1 

0 

0 

1 

1 

0 

0 

0 

0 

0 

.  1 

0 

0 

0 

0 

0 

0 

-71  -  72 

0 

0 

0 

0 

0 

1 

0 

-72 

0 

0 

0 

0 

0 

0 

0 

-2/i 

0 

0 

0 

0 

0 

0 

0 

-/* 

1 

0 

0 

0 

0 

1 

-1 

-fi-fl 

0 

0 

0 

0 

0 

0 

0 

-/* 

1 

0 

0 

0 

k 

0 

-1-ib 

-M-72  . 

Equation  55  can  be  solved  by  standard  techniques  to  yield  all  the  time  evolving  second 
order  moments  of  the  system. 


3  Results  and  Figures 

Figure  1  shows  the  result  of  the  analysis  for  the  case  of  71  =  1.0,  72  =  1.0,  =  0.01,  and 

coupling  spring  constant  k  =  1.0.  Figure  2  shows  the  Monte  Carlo  results  for  the  same  system. 

Agreement  between  the  two  sets  of  results  is  generally  very  good  with  no  major  points 
of  difference.  In  general,  the  analytical  results  show  slightly  deemphasized  local  maxima  and 
minima  as  compared  to  the  Monte  Carlo  results.  This  effect  is  most  pronoimced  in  the  cross 
correlation  curves,  such  as  <  >,  Figures  l(c,e,f  j)  and  2(c,e,fJ)  that  have  a  local  maximum 

at  about  5  seconds.  Comparison  of  Figures  1(c),  and  2(c)  shows  this  effect  most  clearly.  It 
is  also  seen  that,  in  general,  the  steady  state  values  of  the  Monte  Carlo  results  are  somewhat 
lower  than  those  of  the  analytical  results.  This  is  consistent  with  the  results  of  the  previous 
two  papers  [2,  3]. 

Four  studies  on  the  effects  of  varying  the  several  parameters,  while  keeping  all  others 
constant,  are  performed:  1)  71  is  varied  from  0.7  to  10,  2)  72  is  varied  from  0.7  to  10.0,  3) 
k,  the  constant  of  the  coupling  spring,  is  varied  from  0.1  to  10.0,  and  4)  Tc  is  increased  from 
0.001  to  10.0. 

The  results  of  the  first  of  these  studies  are  shown  in  Figures  3  and  4.  The  results  are  as 
would  be  expected  from  [2,  3}.  The  responses  of  the  linear  oscillator  do  not  change  shape 
significantly  with  chan^g  values  of  71;  however,  the  magnitude  does  decrease  with  increasing 
values  of  71.  This  is  a  direct  result  of  the  physics  of  the  problem:  increased  damping  of  an 
oscillator  should  decrease  the  excursions.  The  decrease  in  the  magnitude  of  the  response  of 
the  Duffing  oscillator  is  a  direct  result  of  the  decreased  magnitude  of  the  response  of  the  linear 
oscillator.  The  Duffing  oscillator  is  excited  by  the  displacement  of  the  linear  oscillator;  it  is 
reasonable  that  a  decrease  in  the  magnitude  of  the  displacement  of  the  linear  osdllator  will  be 
reflected  in  a  decrease  in  the  response  of  the  Duffing  osdllator.  The  agreement  between  the 
analytical  and  numerical  data  is  good.  The  analytical  data  show  slightly  deemphasized  local 
maxima  and  minima,  although  the  difference  is  not  marked.  The  largest  difference  is  seen  in 
Figures  3(j)  and  4(j).  Comparison  of  these  two  figure  shows  the  Monte  Carlo  traces  becoming 
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negative  at  five  and  ten  seconds  whUe  the  analytical  traces  merely  reach  local  minima  of  about 
zero. 

The  results  of  the  study  of  the  effect  of  varying  73  is  shown  in  Figures  5  and  6.  These  figures 
show  no  change  in  the  response  of  the  linear  oscillator  with  changing  73.  This  is  because,  while 
the  Duffing  oscillator  is  excited  by  the  displacement  of  the  linear  oscillator,  the  linear  oscillator 
is  independent  of  the  Duffing  oscillator;  i.e.,  neither  nor  03  appear  in  the  equations  of  motion 
of  the  linear  oscillator.  Therefore,  one  would  not  expect  any  variation  of  the  Duffing  oscillator 
to  affect  the  linear  oscillator.  However,  the  effect  on  the  Duffing  oscillator  of  varying  72  is 
quite  pronounced  and  similar  to  the  effect  on  the  linear  oscillator  of  varying  71.  All  response 
curves  of  the  Duffing  osdllator  decrease  in  magnitude  as  73  increases.  This  is  for  the  reasons 
stated  above.  Comparison  of  the  Figures  5  and  6  shows  the  same  similarities  and  contrasts  as 
for  the  7i  study. 

Figures  7  and  8  show  the  results  of  varying  the  spring  constant,  k,  of  the  coupling  spring. 
Again  it  is  seen  that  varying  the  nature  of  the  coupling  has  no  effect  on  the  linear  oscillator: 
there  is  no  change  in  the  moments  of  the  linear  oscillator,  <  Ci  >>  <  Ci*?i  ><  <  rH  >.  The 

moments  of  the  Duffing  oscillator,  <  C3  >t  <  0*72  >>  ^d  <  >,  do  increase  with  increasing 

values  of  k.  This  is  again  because  the  driving  force  of  the  Duffing  oscillator  is  l:(Ci  -  Cz)* 
Therefore,  the  larger  values  of  k  will  lead  to  larger  values  of  the  fordng  function,  and  so  larger 
values  of  the  moments.  The  cross  correlations,  <  qi^3  >  <  Ciqz  >•,  Figures  7(f,h)  and  8(f,h) 
show  a  change  in  the  phase  angle  between  the  two  oscillators.  Looking  at  the  progression  of 
traces  as  ib  is  varied  shows  markedly  different  results  than  any  other  presented  in  this  chapter. 
In  all  other  figures  presented,  as  a  parameter  is  varied,  the  curves  change  monotonically  in 
magnitude,  but  do  not  change  significantly  in  shape.  In  these  figures,  the  trace  corresponding 
to  A  =  10.0,  the  largest  value  of  the  parameter,  has  the  median  value.  Also,  the  traces 
corresponding  to  the  two  lowest  values  of  the  parameter.  A:  =  0.1  and  k  =  0.3,  do  not  show 
any  discernible  local  maxima  in  the  first  2.5  seconds;  they  increase  monotonically  from  zero  to 
thdr  majdma.  The  other  three  traces  show  increasing  local  maxima  with  the  local  maximum 
of  the  k  =  10.0  trace  of  roughly  three  times  as  great  as  the  steady  state  value.  The  decrease 
in  the  magnitude  of  these  traces,  after  each  reaches  its  local  maximum,  represents  a  shift  in 
the  phase  angle  between  the  two  oscillators.  Comparison  of  the  two  methods  of  solution  once 
again  shows  very  good  agreement  with  the  same  qualitative  differences  discussed  above. 

Figures  9  and  10  show  the  effects  on  the  response  of  the  oscillators  of  varying  the  correlation 
time  of  the  noise.  As  in  the  previous  papers,  it  is  seen  that  there  is  negligible  effect  in  varying 
Tc  from  0.001  to  0.01,  and  very  small  effect  in  changing  Tc  to  0.1.  For  this  oscillator,  noise  of 
correlation  time  Tc  <  0.01  can  be  considered  as  white  without  any  loss  of  information.  For 
many  systems,  the  white  noise  approximation  is  acceptable  for  the  case  of  colored  noise  with 
correlation  time  Tc  =  0.1.  As  in  the  other  studies,  the  agreement  between  the  analytical  and 
experimental  results  was  quite  good. 
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e.  <zetel  zeU2>  vs  time 


f.  <etal  zeta2>  vs  time 


Figure  3:  Effect  of  increasing  71  on  the  response  of  coupled  oscillators,  72  =  1.0,  k  =  1.0,  Tc  =  0.01 
calculated  analytically. 
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Figure  3  (Continued):  Effect  of  increasing  71  on  the  response  of  coupled  oscillators,  72 
A;  =  1.0,  Te  =  0.01,  calculated  analytically. 
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i.  <etal  eta2>  vs  time 
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Figure  4  (Continued):  Effect  of  increasing  71  on  the  response  of  coi’pled  oscillators,  72  =  1.0, 
fc  =  1.0,  Tc  =  0.01,  calculated  via  a  Monte  Carlo  experiment. 
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e.  <zetol  zett2>  vs  time 


f.  <etal  zela2>  vs  time 


Figure  5:  Effect  of  increasing  72  on  the  response  of  coupled  oscillators,  71  =  1.0,  k  =  1 
calculated  analytically. 
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Figure  5  (Continued);  Effect  of  increasing  72  on  the  response  of  coupled  oscillators,  71 
k  =  1.0,  Tc  =:  0.01,  calculated  analytically. 


Figiire  6:  Effect  of  increasing  73  on  the  response  of  coupled  oscillators,  71  =  1.0,  k  =  1.0,  Tc  =  0.01, 
calculated  via  a  Monte  Carlo  experiment. 
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Figure  6  (Continued):  Effect  of  increasing  73  on  the  response  of  coupled  osallatois  71 
k  1.0,  Te  =  0.01,  calculated  via  a  Monte  Carlo  exp^ment. 
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g.  <eta2^>  vs  time 
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Figure  7  (Continued):  Effect  of  varying  k  on  the  response  of  coupled  oscillators,  71  =  1.0,  72  =  1.0, 
Tc  —  0.01,  calculated  analytically. 
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i.  <eul  eta2>  vs  time 


j.  ‘aeta2  eta2>  vs  time 
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Figure  8  (Continued):  Effect  of  varying  k  on  the  response  of  coupled  oscillators,  71 
Te  =  0.01,  calculated  via  a  Monte  Carlo  experiment. 


Figure  9:  Effect  of  varying  Tc  on  the  response  of  coupled  osdllators,  71  =  1.0,  72  =  1.0,  k 
calculated  analytically. 
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g.  <eU2  >  vs  time 
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Figure  9  (Continued):  Effect  of  varying  Tc  on  the  response  of  coupled  oscillators,  71  =  1.0,  72  =  1.0, 
k  =  1.0,  calculated  analytically. 
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Figure  10:  EflPect  of  varying  Tc  on  the  response  of  coupled  oscillators,  7i  =  1.0,  72  =  1.0,  k  =  1.0, 
calculated  via  a  Monte  Carlo  experiment. 
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Figure  10  (Continued):  Effect  of  varying  Tc  on  the  response  of  coupled  oscillators,  71  =  1.0, 72  =  1.0, 
=  1.0,  calculated  via  a  Monte  Carlo  experiment. 
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4  Conclusions 


The  overwhelming  similarity  between  the  results  given  by  the  two  methods  implies  that 
this  adaptation  of  the  van  Kampen  expansion  is  an  accurate  tool  for  predicting  the 
statistics  of  the  response  of  a  coupled  Duffing-linear  oscillator  excited  by  colored  noise. 
Although  not  directly  shown  here,  it  seems  reasonable  to  assume  that  this  method  would 
give  good  results  for  an  arbitrary  combination  of  linear  and  Duffing  oscillators.  It  was 
also  seen  that,  depending  on  the  magnitude  of  the  correlation  time  as  compared  to  the 
natural  period  of  the  oscillator,  simplifying  assumptions  can  be  made  that  obviate  the 
need  for  this  adaptation. 
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